mathiumdatatest_2_fft.m

详细介绍了FFT在信号处理中的应用

clear
clc
%----------------------------------------
%(1) read data
%(2) FFT method
%
%----------------------------------------
% case 1: y(1,:), stationary data
% case 2: y(2,:), non-stationary data
% case 3: y(3,:), non-stationary data

%  data format:    
        fid1=fopen ('D:\MathieuDataTest\mathieu\del1\time_del1_run01.dat','r');
           y=fscanf(fid1,'%e  ',[3 inf]); 
        fclose(fid1); 

        Ns=length(y);
   
        y1=y(1,:);
        y2=y(2,:);
        y3=y(3,:);

        fs=1000;  % Sampling frequency
        t=0:(1/fs):((Ns-1)/fs);
     
   figure (1)
        subplot(3,1,1);plot(t,y1);
               ylabel ('y1/ ','fontsize',14,'fontweight','bold')
        subplot(3,1,2);plot(t,y2);
               ylabel ('y2 / ','fontsize',14,'fontweight','bold')
        subplot(3,1,3);plot(t,y3);
               ylabel ('y3/ ','fontsize',14,'fontweight','bold')
               xlabel ('Time/ ','fontsize',14,'fontweight','bold')
               
 pause
   
%------------------sepctral computing program hereon ---------------------------------
% 
figure (2)
          freq=(fs/Ns)*[0:Ns/2-1];
          A_y1=fft(y1);
          A_y1=A_y1(1:Ns/2)/Ns*2;
          subplot(3,1,1);plot(freq,abs(A_y1))
          ylabel('Amplitude spectra of y1');xlabel('Frequency/Hz')
%-------------------------          
nlines=128*2;    % It is important for obtained values and points
nft=2*nlines;
nav=fix(Ns/nft);  %Requires NOVERLAP to be strictly less than the window length.
WINDOW = HANNING(nft);

[P,F] = SPECTRUM(y2,y3,nft,nav,WINDOW,fs);
% %   P = SPECTRUM(X,Y) performs spectral analysis of the two sequences
%     X and Y using the Welch method. SPECTRUM returns the 8 column array
%     P = [Pxx Pyy Pxy Txy Cxy Pxxc Pyyc Pxyc]
%     where
%        Pxx  = X-vector power spectral density
%        Pyy  = Y-vector power spectral density
%        Pxy  = Cross spectral density
%        Txy  = Complex transfer function from X to Y = Pxy./Pxx
%        Cxy  = Coherence function between X and Y = (abs(Pxy).^2)./(Pxx.*Pyy)
%        Pxxc,Pyyc,Pxyc = Confidence range.
%
% % % figure (3)
% % % subplot(4,1,1);plot(F,Pxx,'-k')
% % % subplot(4,1,2);plot(F,Pyy,'-b')
% % % subplot(4,1,3);plot(F,abs(Txy),':r')
% % % subplot(4,1,4);plot(F,angle(Txy),':r')
% % % xlabel ('Frequency/Hz');ylabel ('Spectral density')
Pxx=P(:,1);   %power spectral density of input_data
Pyy=P(:,2);   %power spectral density of output_data
Pxy=P(:,3);   %cross spectral density
Txy=P(:,4);   %Complex transfer function from X to Y = Pxy./Pxx
Cxy=P(:,5);   %Coherence function between X and Y = (abs(Pxy).^2)./(Pxx.*Pyy)
pause
figure (4)   % Y axis should be in log--------------------------
subplot(1,2,1);plot(F,Pxx,'-k');ylabel ('Spectral density of y2');xlabel('Frequency/Hz')
subplot(1,2,2);plot(F,Pyy,'-b');ylabel ('Spectral density of y3');xlabel('Frequency/Hz')
break