image thumbnail

shape analysis of oscillations

by

 

Analysis the of shape of sinusoidal osczillations by interpreting the peaks in the Fourier spectrum.

fft_fit_GUI.m
function [handles,peakhoehe,FT,xmin,xmax,x_GFmin,x_GFmax,xfit,GF,SGFD,fit_data,wid] = fft_fit_GUI(handles,data,schleife,xmin,xmax,x_GFmin,x_GFmax,xfit,GF,SGFD,wid)
% This function contains the main content for calculation.
zeilen = str2num(get(handles.EDIT_zeilen,'String'));
spalte = get(handles.PUM_spalte,'Value');
AF = str2num(get(handles.EDIT_af,'String'));
fs = str2num(get(handles.EDIT_fs,'String'));
anzGF = str2num(get(handles.EDIT_anzGF,'String'));
x0 = ones(zeilen,1);
AHH_eingabe = str2num(get(handles.EDIT_ahh,'String'));
count_sin = handles.count_sin;
   
%%-Auswahl der FFT-Analyse-Methode: 1 cutoff//(2,3,4) Fensterfunktionen
switch get(handles.PUM_rauschen,'Value')
    case 1  %--cutoff------------------------------------------------------
     %% finden des FT-Spektrums mit dem geringsten Rauschen
    for i=1:handles.cutoff
        x=data(i:zeilen,spalte);
        if i==1
            axes(handles.AXES_sin)
            plot(x);   % x plotten
            count_sin = count_sin + 1;
            xlabel(handles.language.plots.xlabel_sin);
            xlim([0 zeilen*get(handles.SLID_sin_x,'Value')]);
        end  
        label = get(handles.PUM_spalte,'String');
        title([handles.language.plots.file,num2str(schleife),': ',char(label(get(handles.PUM_spalte,'Value')))]);
        
        m = length(x);          % Window length
        %n = pow2(nextpow2(m)); % Transform length
        y = fft(x-mean(x),m);   % DFT
        f = (0:m-1)*(fs/m);     % Frequency range
        power = y.*conj(y)/m;   % Power of the DFT
        maximum=max(y(2:length(y)));
        normy=y/maximum;
        
        
        axes(handles.AXES_fft); % AXES_fft als aktive AXES festlegen
        if schleife==1 && i==1  % Bereich der Minimumsbestimmung festlegen
            
            Mitte = fs/2;           % Mitte des Spektrums, entspricht der Nyquistfrequenz
            MitteD = zeilen/2;      % Mitte des Spektrums, in Datenpunkten
            AHH = Mitte/AF;         % Anzahl an hheren Harmonischen
            SAFD = MitteD/AHH;      % Stelle in Datenpunkten, an der die Anregefrequenz liegt
            
            help_x = 0:m-1;
            semilogy(help_x,abs(normy),'b');
            axis([0 (AHH_eingabe+1)*SAFD 1e-7 1]);  
            title(handles.language.plots.titlefft1);            
            xlabel(handles.language.plots.xlabel_sin);
            [x,y] = ginput(2);
            xmin=round(x(1));
            xmax=round(x(2));
            hold on
            semilogy(help_x(xmin:xmax),abs(normy(xmin:xmax)),'g');
            
            
            title(handles.language.plots.titlefft2);
            [xGF,y] = ginput(1);  
           % SGFDg = (xGF/AF)*SAFD;    % Stelle der gesuchten Frequenz in Datenpunkten
                x_GFmin(1)=round(xGF-0.25*SAFD);   % 0.3% entspricht bei 15000 Zeilen 45 Datenpunkten
                x_GFmax(1)=round(xGF+0.25*SAFD);
          %  x_GF(1)=round(x(1)-0.003*zeilen);   % 0.3% entspricht bei 15000 Zeilen 45 Datenpunkten
          %  x_GF(2)=round(x(1)+0.003*zeilen);
            title('');
        end
    
        semilogy(f,abs(normy),'b');
        axis([0 ((AHH_eingabe+1)*AF) 1e-7 1]);
        hold on
        semilogy(f(xmin:xmax),abs(normy(xmin:xmax)),'g');   %farbliche kennzeichnung des Bereiches, wo die
        semilogy(f(x_GFmin(1):x_GFmax(1)),abs(normy(x_GFmin(1):x_GFmax(1))),'r');
        minmean(i,1)=mean(abs(normy(xmin:xmax)));             % Bestimmung des Minimums stattfindet
        minmean(i,2)=mean(power(xmin:xmax));             % Bestimmung des Minimums stattfindet   
        grid on
        
        % ber Anregefrequenz AF, zu untersuchender Frequenz GF Anzal an zeilen z.B. 60000 und fs z.B. 1000
        % wird nun der Bereich des gewnschten peaks berechnet
        
     %   Mitte = fs/2;           % Mitte des Spektrums, entspricht der Nyquistfrequenz
     %   MitteD = zeilen/2;      % Mitte des Spektrums, in Datenpunkten
     %   AHH = Mitte/AF;         % Anzahl an hheren Harmonischen
     %   SAFD = MitteD/AHH;      % Stelle in Datenpunkten, an der die Anregefrequenz liegt
     %   SGFD = (GF/AF)*SAFD;    % Stelle der gesuchten Frequenz in Datenpunkten
        
        
        hold off
        title1 = get(handles.PUM_spalte,'String');
        title(title1(get(handles.PUM_spalte,'Value')));
        xlabel('frequency [Hz]');
        ylabel('peakintensity');
        drawnow
    end
    [val,minmean_idx]=min(minmean(:,2));
    minmean_idx=minmean_idx(1); %selection for multiple minimas
    
    %% Calculating the FFT for a length with minimum mean background

    r=minmean;               % i=minmean % i wurde hier in r getauscht, da i schon verwendet wird
    x=data(minmean_idx:zeilen,spalte);
    m = length(x);          % Window length
    %n = pow2(nextpow2(m));  % Transform length
    y = fft(x,m);           % DFT
    f = (0:m-1)*(fs/m);     % Frequency range
    power = y.*conj(y)/m;   % Power of the DFT
    maximum=max(y(2:length(y)));
    normy=y/maximum;
 %   semilogy(f,abs(normy),'b'); % FFT mit geringstem Rauschen plotten
 %   axis([0 fs/2 1e-7 1]);
        
 
    
    case {2,3,4}  %--Fensterfunktionen-------------------------------------
        x=data(1:zeilen,spalte);
        axes(handles.AXES_sin)
        plot(x);   % x plotten
        count_sin = count_sin + 1;
        xlabel(handles.language.plots.xlabel_sin);
        label = get(handles.PUM_spalte,'String');
        title([handles.language.plots.file,' ',num2str(schleife),': ',char(label(get(handles.PUM_spalte,'Value')))]);
        xlim([0 zeilen*get(handles.SLID_sin_x,'Value')]);
        switch get(handles.PUM_rauschen,'Value')
            case 2
                % Hanning-Window
                x = hann(zeilen).*x;
            case 3
                % Dreiecksfenster
                x = triang(zeilen).*x;
            case 4
                % Blackman-Window
                x = blackman(zeilen).*x;
        end
        m = length(x);  % Window length
        %n = pow2(nextpow2(m)); % Transform length
        y = fft(x,m);           % DFT
        f = (0:m-1)*(fs/m);     % Frequency range
        power = y.*conj(y)/m;   % Power of the DFT
        
        Mitte = fs/2;           % Mitte des Spektrums, entspricht der Nyquistfrequenz
        MitteD = zeilen/2;      % Mitte des Spektrums, in Datenpunkten
        AHH = Mitte/AF;         % Anzahl an hheren Harmonischen
        SAFD = MitteD/AHH;      % Stelle in Datenpunkten, an der die Anregefrequenz liegt
        
        %max(y(max(SAFD-20,1):min(SAFD+20,length(x))))
        maximum=max(y(SAFD-20:SAFD+20));                      % BERARBEITEN!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        %maximum=max(y(25:end-25));                      % BERARBEITEN!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        normy=y/maximum;
       % semilogy(abs(normy))
       % axis([0 zeilen/200 1e-8 100]);
        
        axes(handles.AXES_fft); % AXES_fft als aktive AXES festlegen
        
        if (schleife == 1) || (get(handles.PUM_fit,'Value')==4)
            semilogy(f,abs(normy),'b');
            axis([0 ((AHH_eingabe+1)*AF) 1e-7 1]);  % (AHH_eingabe+1), damit die AHH_eingabe noch angezeigt wird und nicht genau auf dem Rand liegt
        end
        if (schleife == 1) && (get(handles.PUM_fit,'Value')~=4)
            for i=1:anzGF
                title([handles.language.plots.titlefft2,' (',num2str(i),'):']);
                xlabel('frequency [Hz]');
                [x,y] = ginput(1);      
                xx(i) = x;
               % Mitte = fs/2;           % Mitte des Spektrums, entspricht der Nyquistfrequenz
               % MitteD = zeilen/2;      % Mitte des Spektrums, in Datenpunkten
               % AHH = Mitte/AF;         % Anzahl an hheren Harmonischen
               % SAFD = MitteD/AHH;      % Stelle in Datenpunkten, an der die Anregefrequenz liegt
                SGFDg = (xx(i)/AF)*SAFD;    % Stelle der gesuchten Frequenz in Datenpunkten
                x_GFmin(i)=round(SGFDg-0.25*SAFD);   % SAFD entspricht Abstand zwischen 2 Hheren Harmonischen
                x_GFmax(i)=round(SGFDg+0.25*SAFD);
                
                hold on
                semilogy(f(x_GFmin(i):x_GFmax(i)),abs(normy(x_GFmin(i):x_GFmax(i))),'r');
            end
            %title('');
        end
        
        if anzGF > 1
            figure(handles.handle_anzGF)
            subplot(anzGF,2,2)
            semilogy(f,abs(normy),'b');
            axis([0 ((AHH_eingabe+1)*AF) 1e-7 1]);
            for i = 1:anzGF
                hold on
                semilogy(f(x_GFmin(i):x_GFmax(i)),abs(normy(x_GFmin(i):x_GFmax(i))),'r');
            end
            title(['Datei ',num2str(schleife)]);
        end        
        axes(handles.AXES_fft); % AXES_fft als aktive AXES festlegen
   %     hb = figure(); 
    %{    
        semilogy(f,abs(normy),'b'); 
        hold on
        semilogy(f(x_GFmin(1):x_GFmax(1)),abs(normy(x_GFmin(1):x_GFmax(1))),'r');
        axis([0 fs/2 1e-7 1]);
        grid on
    %}    
   %     axes(handles.AXES_fft); % AXES_fft als aktive AXES festlegen
end

    
    if (schleife == 1) && (get(handles.PUM_fit,'Value')~=4)
        for i=1:anzGF
            [val_GF_i,posi_GF] = max(abs(normy(x_GFmin(i):x_GFmax(i))));    % posi_GF ist die Stelle des Maximums aus dem Bereich, in dem GF liegt 
            val_GF(i) = val_GF_i;
            SGFD(i) = posi_GF+x_GFmin(i)-1; % Stelle der gesuchten Frequenz in Datenpunkten%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            Mitte = fs/2;           % Mitte des Spektrums, entspricht der Nyquistfrequenz
            MitteD = zeilen/2;      % Mitte des Spektrums, in Datenpunkten
            AHH = Mitte/AF;         % Anzahl an hheren Harmonischen
            SAFD = MitteD/AHH;      % Stelle in Datenpunkten, an der die Anregefrequenz liegt
            GF(i) = SGFD(i)/SAFD*AF;
        end
    end
    if ((get(handles.PUM_fit,'Value')==2) || (get(handles.PUM_fit,'Value')==3))...
                && schleife == 1
        set([handles.EDIT_pause,handles.SLIDER_pause,handles.TEXT_pause],'Visible','on');
    end

    % y-werte des FFT-spektrums mit dem geringsten Rauschen
    y=abs(normy);
    fr=f';      
    laenge=length(fr);
    halbelaenge=laenge/2;
    FT = [fr(1:halbelaenge),y(1:halbelaenge)];
    
    %% peakhhe bestimmen
    peakhoehe = ones(1,anzGF);
    switch get(handles.PUM_fit,'Value')
        case 1
            %% mit max-Funktion
            for i=1:anzGF
                [peakhoehe_i,posi]= max(y(SGFD(i)-20:SGFD(i)+20));
                peakhoehe(i) = peakhoehe_i;
            %[peakhoehe,posi]=
            %max(y(max(SGFD-20,1):min(SGFD+20,length(x))))  % macht Stress
                    %mit der Fensterfunktion - why ever -.-
                xfit{i}(1) = peakhoehe(i);
                xfit{i}(2) = GF(i); % BERARBEITEN!!
                xfit{i}(3) = -1;  xfit{i}(4) = -1; xfit{i}(5) = -1;% Null-Werte, da keine Daten durch max-Funktion
            end        
            fit_data{1} = 0;
            fit_data{2} = 0;
            fit_data{3} = 0;
                        
        case {2,3}
            %% mit xfit
            xdata = cell(1,anzGF);
            ydata = cell(1,anzGF);            
            ytest = cell(1,anzGF);
            
            % Zu untersuchende Breite um den Peak bestimmen
            if schleife == 1
                %-Peakbreite abschtzen (1. Gesuchte Frequenz)
                r = 1;      % r = Breite/Sigma in Datenpunkten 
                while y(SGFD(1)+r)-y(SGFD(1)+r-1) > y(SGFD(1)+r+1)-y(SGFD(1)+r)
                    r = r+1;
                end
                wid = 20*r; % wid = halbe Fensterbreite in Datenpunkten
            end
       
            for i=1:anzGF
                xdata{i}=f(SGFD(i)-wid:SGFD(i)+wid);          % der variable xdata werden die x-werte +/- wid zugewiesen
                xdata{i}=xdata{i}';
                ydata{i}=abs(normy(SGFD(i)-wid:SGFD(i)+wid)); % der varibale ydata werde die zugehrigen y-werte zugewiesen
                % Startwerte fr xfit:
                if schleife == 1 
                    xfit{i}(1) = val_GF(i);        % Amplitude 
                    xfit{i}(2) = GF(i);            % Erwartungswert
                    %-Breite abschtzen
                    r1 = 0;     %-rechte Seite                 % Breite/Sigma in Datenpunkten 
                    while (y(SGFD(i)+r1+1)-y(SGFD(i)+r1) < 0) && (r1<200)  %y(SGFD(i)+r1)-y(SGFD(i)+r1-1) > y(SGFD(i)+r1+1)-y(SGFD(i)+r1) %y(SGFD(i)+r1+1)-y(SGFD(i)+r1) < 0  % y(SGFD(i)+r)-y(SGFD(i)+r-1) > y(SGFD(i)+r+1)-y(SGFD(i)+r)  r0 = 1!!
                        r1 = r1+1;
                    end  
                    r2 = 0;     %-linke Seite                         % Breite/Sigma in Datenpunkten 
                    while (y(SGFD(i)-r2-1)-y(SGFD(i)-r2) < 0) && (r2<200)  
                        r2 = r2+1;
                    end
                    r = (r1+r2)/6;  %=mean(r1,r2)/3
                    xfit{i}(3) = r*fs/zeilen;      % Breite/sigma
%                     disp(xfit{i}(3))
                    yrausch = [ydata{i}(1:wid+1-4*r);ydata{i}(wid+1+4*r:end)];    % Rauschen ohne die Punkte um den Peak
                    xfit{i}(4) = mean(yrausch);    % Rauschen
                end

                % FITTEN
                % myfunc ist eine funktion mit der gau-fitfunktion die in einem m-file im
                % arbeitsordner gespeichert sein muss
                % 0.0005 ist die amplitude, also die flche unter dem peak
                % 40 ist die frequenz, wo der peak erwartet wird
                % 0.05 ist das sigma, also die peakbreite
                % 0.0002 ist das m, also das level des hintergrundrauschens
                xfit_i = xfit{i};
                xdata_i = xdata{i};
                ydata_i = ydata{i};
            
                [xfit_i] = lsqcurvefit(@mygauss,xfit_i, xdata_i, ydata_i);      % vllt mit 'TolFun' die Toleranzen verndern, aber glaub nicht ntig
                xfit{i}=xfit_i;
                % in [xfit] werden vier werte gespeichert von denen die letzten zwei die
                % peakbreite und das hintergrundrauschen sind. die ersten beiden sind peak-hhe und
                % peak-position
%                 disp(xfit_i(3))
                        
                if anzGF > 1
                    figure(handles.handle_anzGF)
                    subplot(anzGF,2,2*i-1)
                    hold off
                    semilogy(xdata{i},ydata{i});
                end
                if i == 1
                    axes(handles.AXES_fft)
                    hold off
                    semilogy(xdata{1},ydata{1}); 
                    ymax = 0.2;
                end

                % ytest sind die y-werte der Fit-funktion
                ytest{i} = (xfit{i}(1)*exp((-(xdata{i}-xfit{i}(2)).^2)./(2.*xfit{i}(3).^2))) + xfit{i}(4);    % Gleiche Funktion auch in PB_fft_plus_Callback @main_window.m
                xfit{i}(5) = trapz(xdata{1},ytest{i}- xfit{i}(4));   % Flche unter dem Graphen
                if anzGF > 1
                    figure(handles.handle_anzGF)
                    hold on
                    subplot(anzGF,2,2*i-1)
                    semilogy(xdata{i},ytest{i},'black');
                    grid on
                    xlabel('frequency [Hz]');
                    ylabel('paekintensity');
                    xlim([GF(i)-wid/2*(AF/SAFD) GF(i)+wid/2*(AF/SAFD)]);        % anpassen?
                    maxytest = max(ytest{i});                    
                    while (maxytest > ymax/10)
                        ymax = ymax * 10;
                    end
                    ylim([0.000001 ymax]);
                end  
                if i == 1
                    axes(handles.AXES_fft)
                    hold on
                    semilogy(xdata{1},ytest{1},'black');
                    if schleife == 1    % Kontrolle, ob der richtige Punkt herausgefunden wurde beim Anklicken
                        plot(f(SGFD(i)),val_GF(i),'ro') 
                    end
                    title1 = get(handles.PUM_spalte,'String');
                    title([char(title1(get(handles.PUM_spalte,'Value'))),...
                        ': peakfitting - ',handles.language.plots.file,' ',num2str(schleife)]);
                    grid on
                    xlabel('frequency [Hz]');
                    ylabel('paekintensity');
                    handles.xlim_fit = [GF(1)-wid/2*(AF/SAFD) GF(1)+wid/2*(AF/SAFD)];
                    xlim([GF(1)-wid/2*(AF/SAFD) GF(1)+wid/2*(AF/SAFD)]);
                    maxytest = max(ytest{1}); 
                    while (maxytest > ymax/10)
                        ymax = ymax * 10;
                    end
                    ylim([0.000001 ymax]);
                end
                hold off

                % maximum in xfit bestimmen
                %peakhoehe ist die hoehe aus xfit
                %posi ist die position, nur noetig, damit funktion mit max funktioniert
                switch get(handles.PUM_fit,'Value')
                    case 2
                        [peakhoehe_i,posi]=max(ytest{i});
                        peakhoehe(i) = peakhoehe_i;
                    case 3
                        [peakhoehe_i,posi]= max(y(SGFD(i)-20:SGFD(i)+20));
                        peakhoehe(i) = peakhoehe_i;                        
                end
            end
            
            fit_data{1} = xdata{1};        % cell-Array in dem die Werte der fit-Kurve gespeichert werden
            fit_data{2} = ydata{1};        % um sie spter durch +/- durchblttern zu knnen
            fit_data{3} = ytest{1};        % ACHTUNG: Nur erste Frequenz wird gespeichert! 
            
        case 4
            %%-Integral----------------------------------------------------
            %semilogy(f,y,'b')
            %axis([0 fs/2 1e-7 1]);
            area = trapz(f,y)
            
            %rauschen = mean(y(SGFD(1):round(end/2)));
            %hold on
            %semilogy(f,rauschen,'k')
            
            %-Nullwerte fr bergabeparamter setzen------------------------
            for i=1:anzGF
                xfit{i}(1) = area;
                xfit{i}(2) = -1; 
                xfit{i}(3) = -1;  
                xfit{i}(4) = -1;
            end        
            fit_data{1} = 0;
            fit_data{2} = 0;
            fit_data{3} = 0;
            
    end
    handles.count_sin = count_sin;
end

Contact us