matlab谐波分析程序

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=1000;
order=7;
%**********************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。