matlab諧波分析程序(共5頁)

1、clcclear all;format long;Ns=1000;order=13; %*read the position and flux density*fid=fopen(B.dat,r); %open the original filefidnew = fopen(b1.dat,w); %write the new filewhile feof(fid)=0 tline = fgetl(fid); %tline? if ischar(tline), break, end temp=abs(tline); Nlength=length(tline); isemptyline=0; %?

2、 if Nlength=0 isemptyline=1; end allspace=0; %? isspace=0; for i=1:Nlength T=temp(i); if T=32 isspace=isspace+1; end if isspace=Nlength allspace=1; break end end findalpha=0; %? for j=1:Nlength T=temp(j); if (T=65)&(T=90)|(T=97)&(T=122) findalpha=1; break; end end if (findalpha)&(allspace)&(isemptyl

3、ine=0) %? fprintf(fidnew,tline);fprintf(fidnew,n); end end fclose(fid);fclose(fidnew); fid1=fopen(b1.dat,r);flux_position =fscanf(fid1,%f,2,Ns);fclose(fid1);%*read file finish*flux_position=flux_position;pos1=flux_position(:,1);pos_delta=pos1(2);pos_length=length(pos1);pos_last=pos1(pos_length);for

4、i=1:1:pos_length %copy and get another part of position pos2(i)=pos_last+i*pos_delta;endpos1=pos1;flux1=flux_position(:,2);flux2=-flux_position(:,2);pos=pos1,pos2;%combine and get all part of positionflux1=flux1;flux2=flux2;flux=flux1,flux2;%combine and get all part of flux density valuefigure;plot(

5、pos1,flux1,r);%plot origional waveformhold on;grid on;fft1=fft(flux,Ns);j=0;amp_har=zeros(1,(order+1)/2); for m=1:2:order j=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

6、_back; fft1=ifft(fft1,Ns); fft1=real(fft1); plot(pos1,fft1); hold on;end k=(1:2:order);figure;bar(k,amp_har);grid on;peak_b=max(fft1)rms_b=0.707*peak_bclcclear all;format long; Ns=1000;order=7; %*read the position and flux density*fid=fopen(B.dat,r); %open the original filefidnew = fopen(b1.dat,w);

7、%write the new filewhile feof(fid)=0 tline = fgetl(fid); %tline? if ischar(tline), break, end temp=abs(tline); Nlength=length(tline); isemptyline=0; %? if Nlength=0 isemptyline=1; end allspace=0; %? isspace=0; for i=1:Nlength T=temp(i); if T=32 isspace=isspace+1; end if isspace=Nlength allspace=1; b

8、reak end end findalpha=0; %? for j=1:Nlength T=temp(j); if (T=65)&(T=90)|(T=97)&(T=122) findalpha=1; break; end end if (findalpha)&(allspace)&(isemptyline=0) %? fprintf(fidnew,tline);fprintf(fidnew,n); end end fclose(fid);fclose(fidnew); fid1=fopen(b1.dat,r);flux_position =fscanf(fid1,%f,2,Ns);fclos

9、e(fid1);%*read file finish*flux_position=flux_position;pos1=flux_position(:,1);pos_delta=pos1(2);pos_length=length(pos1);pos_last=pos1(pos_length);for i=1:1:pos_length %copy and get another part of position pos2(i)=pos_last+i*pos_delta;endpos1=pos1;flux1=flux_position(:,2);flux2=-flux_position(:,2);

10、pos=pos1,pos2;%combine and get all part of positionflux1=flux1;flux2=flux2;flux=flux1,flux2;%combine and get all part of flux density valuefigure;plot(pos1,flux1,r);%plot origional waveformhold on;grid on;fft1=fft(flux,Ns);j=0;amp_har=zeros(1,(order+1)/2); for m=1:2:order j=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