fourierdec.m

关于电机的.m程序

%Fourierdec.m ==> Fourier series decomposition of periodic functions of identical periods 
% specify the column vector t and array f representing the signals over one complete period
% specify the number of harmonics to be calculated by h
% for an illustration, run fourier0.m
clc,close all
disp(date)
disp('fundamental frequency (Hz)')
t = t - t(1); % shift time array such that it starts with 0
L=length(t);
T=t(L);
fo=round(1/T);
disp(fo)
nu=0;
%%%
for w=1:size(f,2);
   nu=nu+1;
   disp(['   Function No. ' num2str(nu)]);
   disp('_________________________________')
   disp('   order of | amplitude |  phase')
   disp('   harmonic |           |  [deg]')
   disp('____________|___________|________')
  
   g=f(:,w);
%%% Complex Fourier coefficients ck
       co=1/T*trapz(t,g)+eps;             
       fprintf('%8.0f%13.4f\n',0,co)
       C=[];
       rms2 = co^2;
     for k=[1:h]
        q=g.*exp(-j*2*pi/T*t*k);  
        ck= 2/T*abs(trapz(t,q));
        pk=angle(trapz(t,q))*180/pi;
        C=[C ck];
        fprintf('%8.0f%13.4f%11.2f\n',k,ck,pk)
        rms2 = rms2+ck^2/2;
     end
     % RMS and THD
     rms = sqrt(rms2);
     THD = sqrt(rms2-C(1)^2/2)/(C(1)/sqrt(2));
     mean = co;
     disp(['RMS:  ' num2str(rms)]);
     disp(['THD:  ' num2str(THD)]);
     disp(['Mean: ' num2str(mean)]);
     
     disp('================================')
    %%%  Plots 
   order=((1:h+1)-1)'; mag=[co C]';    
  figure('Position',[225 70 760 555],'Name','Fourier analysis');
  subplot(2,1,1),plot([0 0],[0 g(1)],'r',t,g,'r',[T T],[g(L) 0],'r',[T 0],[0 0],'k','linewidth',2)
   xlabel('TIME [s]') ;ylabel('SIGNAL');
   title(['PERIODIC SIGNAL (RMS=' num2str(rms) ' THD=' num2str(THD) ')']),grid
   subplot(2,1,2),stem(order,mag,'-ob')
   xlabel('HARMONIC ORDER'); ylabel('AMPLITUDE');
   title('FOURIER COMPONENTS'),grid
end