%Einlesen der Messdaten und die Detektion von Extremstellen

clear all
clc

%#########################################################################
%########  START: Einlesen der 10 Hz Daten aus der einer Stichprobe ######
%#########################################################################
disp('START: Einlesen der 10 Hz Daten aus der einer Stichprobe')
tic;

DELIMITER = '\t';
HEADERLINES = 0;


fid = fopen('.....','r');
data=textscan(fid,'%n ...','delimiter','\t','Headerlines',0); 
fclose(fid);


 ZEIT=data{:,1};

 vh_wind=data{:,4};

 
 Messwerte=[ZEIT vh_wind];
 
ZEIT=Messwerte(:,1);
vh_wind=Messwerte(:,2);


toc;
%#########################################################################
%########  ENDE: Einlesen der 10 Hz Daten aus der einer Stichprobe '######
%#########################################################################


d = sign(diff(vh_wind));                              % Differenz zum VorgŠngerelement bestimmen
f=find(d ~= [d(1); d(1: end - 1)]);                     % Indizes der lokalen Extremwerte berechnen


%#########################################################################
%########  START: grafische Darstellung der Extremwertpositionen #########
%#########################################################################

disp('START: grafische Darstellung der Extremwertpositionen')
tic;

%Ändern der Positionseinträge von vh_wind, wo Extremwerte vorhanden sind. Diese
%werden ersetzt und später durch ein Balkendiagramm ausgegeben. 

M=vh_wind;
z=size(M);
n=z(1,1);
indikator=-1.0000001;   %(Stellgröße) Extremwerte werden durch die Indikatorzahl ersetzt, dies sollte jedoch eine Zahl sein, die im Messsignal nicht erscheinen kann!!!
M(f)=indikator;

INDIKATOR=[];
toc;


disp('Vektor der Länge von vh_wind wird erstellt und die Extremwerte an den jeweiligen Positionen eingetragen')
tic;
for i=1:1:n
    x=M(i);
    if(x==indikator)
        z=2;           %(Stellgröße) wie hoch sollte der Ausschlag angezeigt werden 
    else
        z=0;
    end
    INDIKATOR=[INDIKATOR;z];
end


toc;
%#########################################################################
%########  ENDE: grafische Darstellung der Extremwertpositionen #########
%#########################################################################

%#########################################################################
%###  START: Untersuchung der linken Flanke einer Extremwertpositionen ###
%#########################################################################

disp('START: Untersuchung der linken Flanke einer Extremwertpositionen')
tic;


alpha=300;                     %Beginnt beim n-ten Extremwert, wenn SP am Anfang schon Extremwert aufweist, kann Mittelwert long nicht gebildet werden!
z=size(f);
n=z(1,1);

Matrix=[];
A=[];

Intervall=20;                                   % (Stellgröße) MittelwertIntervall gibt an, wie viele Stellen vom Extremwert aus betrachtet für die Mittelwertbildung verwendet werden sollen
MittelwertIntervall=Intervall-1;                % Später sollte das Intervall ca. 2-3 Sekunden betragen

for i=alpha:1:n                                     %Beginnt beim ersten Extremwert
    for j=MittelwertIntervall:-1:0             % nimmt mir hier 2 Werte vor dem Extremwert und lässt die Schleife rückwärts laufen, hier insgesamt 3 Werte
    extremwert2=vh_wind(f(i)-j); 
    A=[A; extremwert2]; 
    end

    Matrix=[A];
end


% Bilden des Mittelwertes an der Peakstelle

PEAKmean=[];

%Intervall=3;                                   % (Stellgröße) MittelwertIntervall gibt an, wie viele Stellen vom Extremwert aus betrachtet für die Mittelwertbildung verwendet werden sollen
MittelwertIntervall=Intervall-1;                % Später sollte das Intervall ca. 2-3 Sekunden betragen

for i=1:MittelwertIntervall+1:length(Matrix)
    x=Matrix(i:i+MittelwertIntervall);        %Mittelwertbildung über 3 Werte
    y=mean(x);
    PEAKmean=[PEAKmean; y];

end

%##########################################################################
%Bilden des Mittelwertes an der Stelle x_____________I__PEAK


z=size(f);
n=z(1,1);

Matrix_long=[];
A_long=[];

t=600;                        % (Stellgröße) gibt das Zeitintervall an, was zwischen x______________I liegt (600 Werte entspricht 1 min bei 10 Hz Daten)
for i=alpha:1:n
    for j=MittelwertIntervall+t:-1:MittelwertIntervall+1             % nimmt mir hier 2 Werte vor dem Extremwert und lässt die Schleife rückwärts laufen, hier insgesamt 3 Werte
    extremwert_long=vh_wind(f(i)-j); 
    A_long=[A_long; extremwert_long]; 
    end

    Matrix_long=[A_long];
end
   

% Bilden des Mittelwertes im Intervall x______________I

PEAKmean_long=[];


Intervall_2=t;                                   % MittelwertIntervall gibt an, wie viele Stellen vom Extremwert aus betrachtet für die Mittelwertbildung verwendet werden sollen
MittelwertIntervall=Intervall_2-1;                % Später sollte das Intervall ca. 2-3 Sekunden betragen

for i=1:MittelwertIntervall+1:length(Matrix_long)
    x=Matrix_long(i:i+MittelwertIntervall);        %Mittelwertbildung über 3 Werte
    y=mean(x);
    PEAKmean_long=[PEAKmean_long; y];
end
    


%%Differenzmatrix=[PEAKmean PEAKmean_long]

differenz=abs(PEAKmean-PEAKmean_long);           % Es wird die absolute Abweichung verwendet, damit auch Flauten und nicht nur Böen detektiert werden können

z=size(differenz);
n=z(1,1);
Schranke=30;        % (Stellgröße) Später bei den Windmessdaten wird die minimale Schranke bei ca. 4,5 m/s liegen
BOENINDI=[];
for i=1:1:n
    x=differenz(i);
    if(x>Schranke)
        y=1;
    else
        y=0;
    end
    BOENINDI=[BOENINDI;y];
end



toc;


%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


%Ändern der Positionseinträge von vh_wind, wo Extremwerte vorhanden sind. Diese
%werden durch die Differenzbeträge der Mittelwerte ersetzt.

disp('Extremwerte werden durch eine Indikatorzahl belegt, die später gefiltert werden kann')
tic;

D=vh_wind;
z=size(D);
n=z(1,1);

for i=1:1:length(differenz)
    indikator=differenz(i);   %(Stellgröße) Extremwerte werden durch die Indikatorzahl ersetzt, dies sollte jedoch eine Zahl sein, die im Messsignal nicht erscheinen kann!!!
    D(f(i))=indikator;         %erstetzt die Extremwertstelle f(i) im kopierten Windvektor durch die Differenz der Mittelwerte (PEAK & long)
end


%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++