fft_bn_3(傅立叶分解m文件）.m

傅立叶分解M文件

n=1; %谐波次数
n2=2
n3=3
format long;
tp=T(:,1);                               %导入的文件名应改为T
yt=T(:,2);
%tp=Bn(:,1);                           %导入的文件名应改为T
%yt=Bn(:,2);
f=50;                   %计算气隙磁没用
%f=35;
T1=1/f;
m=(tp(length(tp))-tp(1))/T1;           %波形包括m个周期
w=2*pi*f;                              %角频率
yt_Bav=abs(yt);                        %将径向磁密取绝对值，为计算平均值用
Bav(1)=(tp(2)-tp(1))*(yt_Bav(1)+yt_Bav(2))/2;
a0(1)=(tp(2)-tp(1))*(yt(1)+yt(2))/2;
an(1)=(tp(2)-tp(1))*(yt(1)*cos(n*w*tp(1))+yt(2)*cos(n*w*tp(2)))/2;  %%%%%开始计算不同频率%%%
bn(1)=(tp(2)-tp(1))*(yt(1)*sin(n*w*tp(1))+yt(2)*sin(n*w*tp(2)))/2;  %%%%%%%%%%开始计算不同频率%%%
for i=1:length(tp)-1
a0(i+1)=a0(i)+(tp(i+1)-tp(i))*(yt(i)+yt(i+1))/2;
an(i+1)=an(i)+(tp(i+1)-tp(i))*(yt(i)*cos(n*w*tp(i))+yt(i+1)*cos(n*w*tp(i+1)))/2;
bn(i+1)=bn(i)+(tp(i+1)-tp(i))*(yt(i)*sin(n*w*tp(i))+yt(i+1)*sin(n*w*tp(i+1)))/2;
Bav(i+1)=Bav(i)+(tp(i+1)-tp(i))*(yt_Bav(i)+yt_Bav(i+1))/2;
end
a0=a0(length(tp))*2/T1/m;
an=an(length(tp))*2/T1/m;
bn=bn(length(tp))*2/T1/m;
An=sqrt(an^2+bn^2)                                 %基波幅值
phA=-atan(bn/an)*180/pi                            %初始相位
%t=0:0.001:length(tp);
Ft1=-An*cos(n*w*tp+phA*pi/180);

an(1)=(tp(2)-tp(1))*(yt(1)*cos(n2*w*tp(1))+yt(2)*cos(n2*w*tp(2)))/2;
bn(1)=(tp(2)-tp(1))*(yt(1)*sin(n2*w*tp(1))+yt(2)*sin(n2*w*tp(2)))/2;
for i=1:length(tp)-1
a0(i+1)=a0(i)+(tp(i+1)-tp(i))*(yt(i)+yt(i+1))/2;
an(i+1)=an(i)+(tp(i+1)-tp(i))*(yt(i)*cos(n2*w*tp(i))+yt(i+1)*cos(n2*w*tp(i+1)))/2;
bn(i+1)=bn(i)+(tp(i+1)-tp(i))*(yt(i)*sin(n2*w*tp(i))+yt(i+1)*sin(n2*w*tp(i+1)))/2;
Bav(i+1)=Bav(i)+(tp(i+1)-tp(i))*(yt_Bav(i)+yt_Bav(i+1))/2;
end
a0=a0(length(tp))*2/T1/m;
an=an(length(tp))*2/T1/m;
bn=bn(length(tp))*2/T1/m;
An=sqrt(an^2+bn^2)                                 %基波幅值
phA=-atan(bn/an)*180/pi                            %初始相位
%t=0:0.001:length(tp);
Ft2=-An*cos(n2*w*tp+phA*pi/180);

an(1)=(tp(2)-tp(1))*(yt(1)*cos(n3*w*tp(1))+yt(2)*cos(n3*w*tp(2)))/2;
bn(1)=(tp(2)-tp(1))*(yt(1)*sin(n3*w*tp(1))+yt(2)*sin(n3*w*tp(2)))/2;
for i=1:length(tp)-1
a0(i+1)=a0(i)+(tp(i+1)-tp(i))*(yt(i)+yt(i+1))/2;
an(i+1)=an(i)+(tp(i+1)-tp(i))*(yt(i)*cos(n3*w*tp(i))+yt(i+1)*cos(n3*w*tp(i+1)))/2;
bn(i+1)=bn(i)+(tp(i+1)-tp(i))*(yt(i)*sin(n3*w*tp(i))+yt(i+1)*sin(n3*w*tp(i+1)))/2;
Bav(i+1)=Bav(i)+(tp(i+1)-tp(i))*(yt_Bav(i)+yt_Bav(i+1))/2;
end
a0=a0(length(tp))*2/T1/m;
an=an(length(tp))*2/T1/m;
bn=bn(length(tp))*2/T1/m;
An=sqrt(an^2+bn^2)                                 %基波幅值
phA=-atan(bn/an)*180/pi                            %初始相位
%t=0:0.001:length(tp);
Ft3=-An*cos(n3*w*tp+phA*pi/180);

plot(tp,yt,'r',tp,Ft1,'k',tp,Ft2,'b',tp,Ft3);
grid;
xlabel('x/mm');
ylabel('B/T');

B_max=max(yt)
B_min=min(yt);
B_max_abs=(B_max+abs(B_min))/2;                     %取一个极距内气隙磁密最大值
%B_av=Bav(length(tp))/(Bn(length(tp),1)-Bn(1,1))      %实际气隙磁密的平均值
B_av=Bav(length(tp))/(T(length(tp),1)-T(1,1))      %实际气隙磁密的平均值
a_i=B_av/An                                        %计算极弧系数