matlab谐波分析程序.doc

clc clear all; format long; Ns=1000; order=13; %**********************read the position and flux density************************ fid=fopen(B.dat,r); %open the original file fidnew = fopen(b1.dat,w); %write the new file while 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; 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)(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]); 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 i=1:1:pos_length %copy and get another part of position pos2(i)=pos_last+i*pos_delta; end pos1=pos1; flux1=flux_position(:,2); flux2=-flux_position(:,2); pos=[pos1,pos2];%combine and get all part of position flux1=flux1; flux2=flux2; flux=[flux1,flux2];%combine and get all part of flux density value figure; plot(pos1,flux1,r);%plot origional waveform hold 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_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_b clc clear all; format long; Ns=