matlab谐波分析程序
clc
clear all;
format long;
Ns=1000;
order=13;
!
%**********************read the position and flux density************************ fid=fopen('','r'); %open the original file
fidnew = fopen('','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('','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);
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=*peak_b
;
clc
—
clear all;
format long;
Ns=1000;
order=7;
%**********************read the position and flux density************************ fid=fopen('','r'); %open the original file
fidnew = fopen('','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('','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=*peak_b