matlab谐波分析程序

clcclearall;formatlong;Ns=1000;order=13;%**********************readthepositionandfluxdensity************************fid=fopen('B.dat','r');%opentheoriginalfilefidnew=fopen('b1.dat','w');%writethenewfilewhilefeof(fid)==0tline=fgetl(fid);%tline?if~ischar(tline),break,endtemp=abs(tline);Nlength=length(tline);isemptyline=0;%????????????ifNlength==0isemptyline=1;endallspace=0;%????????????isspace=0;fori=1:NlengthT=temp(i);ifT==32isspace=isspace+1;endifisspace==Nlengthallspace=1;breakendendfindalpha=0;%?????????????forj=1:NlengthT=temp(j);if((T=65)&(T=90))|((T=97)&(T=122))findalpha=1;break;endendif(~findalpha)&(~allspace)&(isemptyline==0)%????????????????????????fprintf(fidnew,tline);fprintf(fidnew,'\n');endendfclose(fid);fclose(fidnew);fid1=fopen('b1.dat','r');flux_position=fscanf(fid1,'%f',[2,Ns]);fclose(fid1);%********************************readfilefinish*****************************************flux_position=flux_position';pos1=flux_position(:,1);pos_delta=pos1(2);pos_length=length(pos1);pos_last=pos1(pos_length);fori=1:1:pos_length%copyandgetanotherpartofpositionpos2(i)=pos_last+i*pos_delta;endpos1=pos1';flux1=flux_position(:,2);flux2=-flux_position(:,2);pos=[pos1,pos2];%combineandgetallpartofpositionflux1=flux1';flux2=flux2';flux=[flux1,flux2];%combineandgetallpartoffluxdensityvaluefigure;plot(pos1,flux1,'r');%plotorigionalwaveformholdon;gridon;fft1=fft(flux,Ns);j=0;amp_har=zeros(1,(order+1)/2);form=1:2:orderj=j+1;fft1=fft(flux,Ns);fund_ele_front=fft1(m+1);fund_ele_back=fft1(Ns+1-m);amp_har(j)=(abs(fund_ele_front))/Ns*2;fft1=0*fft1;fft1(m+1)=fund_ele_front;fft1(Ns+1-m)=fund_ele_back;fft1=ifft(fft1,Ns);fft1=real(fft1);plot(pos1,fft1);holdon;endk=(1:2:order);figure;bar(k,amp_har);gridon;peak_b=max(fft1)rms_b=0.707*peak_bclcclearall;formatlong;Ns=1000;order=7;%**********************readthepositionandfluxdensity************************fid=fopen('B.dat','r');%opentheoriginalfilefidnew=fopen('b1.dat','w');%writethenewfilewhilefeof(fid)==0tline=fgetl(fid);%tline?if~ischar(tline),break,endtemp=abs(tline);Nlength=length(tline);isemptyline=0;%????????????ifNlength==0isemptyline=1;endallspace=0;%????????????isspace=0;fori=1:NlengthT=temp(i);ifT==32isspace=isspace+1;endifisspace==Nlengthallspace=1;breakendendfindalpha=0;%?????????????forj=1:NlengthT=temp(j);if((T=65)&(T=90))|((T=97)&(T=122))findalpha=1;break;endendif(~findalpha)&(~allspace)&(isemptyline==0)%????????????????????????fprintf(fidnew,tline);fprintf(fidnew,'\n');endendfclose(fid);fclose(fidnew);fid1=fopen('b1.dat','r');flux_position=fscanf(fid1,'%f',[2,Ns]);fclose(fid1);%********************************readfilefinish*****************************************flux_position=flux_position';pos1=flux_position(:,1);pos_delta=pos1(2);pos_length=length(pos1);pos_last=pos1(pos_length);fori=1:1:pos_length%copyandgetanotherpartofpositionpos2(i)=pos_last+i*pos_delta;endpos1=pos1';flux1=flux_position(:,2);flux2=-flux_position(:,2);pos=[pos1,pos2];%combineandgetallpartofpositionflux1=flux1';flux2=flux2';flux=[flux1,flux2];%combineandgetallpartoffluxdensityvaluefigure;plot(pos1,flux1,'r');%plotorigionalwaveformholdon;gridon;fft1=fft(flux,Ns);j=0;amp_har=zeros(1,(order+1)/2);form=1:2:orderj=j+1;fft1=fft(flux,Ns);fund_ele_front=fft1(m+1);fund_ele_back=fft1(Ns+1-m);amp_har(j)=(abs(fund_ele_front))/Ns*2;fft1=0*fft1;fft1(m+1)=fund_ele_front;fft1(Ns+1-m)=fund_ele_back;fft1=ifft(fft1,Ns);fft1=real(fft1);plot(pos1,fft1);holdon;endk=(1:2:order);figure;bar(k,amp_har);gridon;%peak_b=max(fft1)%rms_b=0.707*peak_b