t=A(J+1);
A(J+1)=A(I+1);
A(I+1)=t;
end
K=N/2;
whileK<=J
J=J-K;
K=K/2;
end
J=J+K;
I=I+1;
end
whileL<=M
LE=2^L;R=1;
LE1=LE/2;
W=1;W1=exp(-j*pi/LE1);
whileR<=LE1
P=R;
whileP<=N
Q=P+LE1;T=A(Q)*W;
A(Q)=A(P)-T;
A(P)=A(P)+T;
P=P+LE;
end
W=W*W1;
R=R+1;
end
L=L+1;
end
end
附录二主程序
closeall;clearall
s_f=15;s_xy=15;s_t=20;
s_mar=10;w_f=2.5;w_axi=2;
%信号产生
T0=8e-2;Ts=1/800;
N=T0/Ts;
t=0:
Ts:
(N-1)*Ts;
f1=50;f2=100;f3=200;fp=80;fs=120;
Sig=4*sin(f1*2*pi*t)+4*sin(f2*2*pi*t)+4*sin(f3*2*pi*t);
s_need=4*sin(f1*2*pi*t);
s_1=4*sin(f2*2*pi*t);
s_2=4*sin(f3*2*pi*t);
P_qian=fftlzr(Sig);
figure
(1)
h1=stem(2*pi/(N*Ts).*(0:
N-1),abs(P_qian));gridon
hx1=xlabel('模拟角频率/rad/s');
ht1=title('滤波前幅度频谱');
ha1=gca;
axistight
%滤波器设计
Fs=1/Ts;
wp=2/Ts*tan(fp*2*pi*Ts/2);
ws=2/Ts*tan(fs*2*pi*Ts/2);
rp=1;rs=15;
[Nl,Wn]=buttord(wp,ws,rp,rs,'s');%求阶数及3dB截止频率
[z,p,k]=buttap(Nl);%求极零点增益
[b,a]=zp2tf(z,p,k);%系统传递函数
[b2,a2]=lp2lp(b,a,Wn);%转换为所需低通滤波器
[HS,WS]=freqs(b2,a2,linspace(0,2000,512));%求0到2000rad/s的频率响应
[hs,wss]=freqs(b2,a2,[wp,ws,[f1,f2,f3].*2*pi]);%求通带阻带边缘角频率及原始信号频率响应
%画模拟滤波器频谱并标注
figure(6)
subplot(211)
hsf=plot(WS,20*log10(abs(HS)));
holdon
stem(wss,20*log10(abs(hs)),'filled','r:
','linewidth',w_f)
text(wss
(1)-100,20*log10(abs(hs
(1)))-10,...
{'\Omega_p=160\pirad/s';[num2str(20*log10(abs(hs
(1)))),'dB']});
text(wss
(2)-65,20*log10(abs(hs
(2)))-12,...
{'\Omega_s=240\pirad/s';[num2str(20*log10(abs(hs
(2)))),'dB']});
text(wss(3)-100,20*log10(abs(hs(3)))-10,...
{'\Omega_1=100\pirad/s';[num2str(20*log10(abs(hs(3)))),'dB']});
text(wss(4)-65,20*log10(abs(hs(4)))-12,...
{'\Omega_2=200\pirad/s';[num2str(20*log10(abs(hs(4)))),'dB']});
text(wss(5)-70,20*log10(abs(hs(5)))-10,...
{'\Omega_3=400\pirad/s';[num2str(20*log10(abs(hs(5)))),'dB']});
hxsf=xlabel('模拟角频率/rad/s');hysf=ylabel('幅值/dB');
htsf=title('幅频特性');hasf=gca;gridon
subplot(212)
hsx=plot(WS,angle(HS));
holdon
stem(wss,angle(hs),'filled','r:
','linewidth',w_f)
text(wss
(1),angle(hs
(1))+0.7,...
{'\Omega_p=160\pirad/s';[num2str(angle(hs
(1))),'rad']});
text(wss
(2),angle(hs
(2))+0.7,...
{'\Omega_s=240\pirad/s';[num2str(angle(hs
(2))),'rad']});
text(wss(3)-150,angle(hs(3))-0.8,...
{'\Omega_1=100\pirad/s';[num2str(angle(hs(3))),'rad']});
text(wss(4)-20,angle(hs(4))+0.8,...
{'\Omega_2=200\pirad/s';[num2str(angle(hs(4))),'rad']});
text(wss(5)-50,angle(hs(5))-1.3,...
{'\Omega_3=400\pirad/s';[num2str(angle(hs(5))),'rad']});
hxsx=xlabel('模拟角频率/rad/s');hysx=ylabel('相值/rad');
htsx=title('相频特性');hasx=gca;gridon
[b3,a3]=bilinear(b2,a2,Fs);%双线性变换到数字滤波器
[H,W]=freqz(b3,a3,linspace(0,0.5*pi,512));%求数字滤波器频率响应
[hd,wd]=freqz(b3,a3,[wp,ws,[f1,f2,f3].*2*pi]*Ts);%求通带阻带边缘数字角频率及原始信号数字频率响应
wd=wd/pi;
%画数字滤波器频谱并标注
figure
(2)
subplot(211)
h2=plot(W/pi,20*log10(abs(H)));
holdon
stem(wd,20*log10(abs(hd)),'filled','r:
','linewidth',w_f)
text(wd
(1)-0.0035,20*log10(abs(hd
(1)))-12,...
{['\omega_p=',num2str(wd
(1)),'\pirad'];[num2str(20*log10(abs(hd
(1)))),'dB']});
text(wd
(2)-0.003,20*log10(abs(hd
(2)))-12,...
{['\omega_s=',num2str(wd
(2)),'\pirad'];[num2str(20*log10(abs(hd
(2)))),'dB']});
text(wd(3)-0.0035,20*log10(abs(hd(3)))-12,...
{['\omega_1=',num2str(wd(3)),'\pirad'];[num2str(20*log10(abs(hd(3)))),'dB']});
text(wd(4)-0.003,20*log10(abs(hd(4)))-12,...
{['\omega_2=',num2str(wd(4)),'\pirad'];[num2str(20*log10(abs(hd(4)))),'dB']});
text(wd(5)-0.003,20*log10(abs(hd(5)))-12,...
{['\omega_3=',num2str(wd(5)),'\pirad'];[num2str(20*log10(abs(hd(5)))),'dB']});
hx2=xlabel('归一化数字角频率/\pirad');
hy1=ylabel('幅值/dB');
ht2=title('幅频特性');
ha2=gca;
gridon
subplot(212)
h3=plot(W/pi,angle(H));grid;title('相频特性');
holdon
stem(wd,angle(hd),'filled','r:
','linewidth',w_f)
text(wd
(1),angle(hd
(1))+1,...
{['\omega_p=',num2str(wd
(1)),'\pirad'];[num2str(angle(hd
(1))),'\pirad']});
text(wd
(2)-0.002,angle(hd
(2))+1,...
{['\omega_p=',num2str(wd
(2)),'\pirad'];[num2str(angle(hd
(2))),'\pirad']});
text(wd(3)-0.003,angle(hd(3))-1,...
{['\omega_p=',num2str(wd(3)),'\pirad'];[num2str(angle(hd(3))),'\pirad']});
text(wd(4)-0.003,angle(hd(4))-2,...
{['\omega_p=',num2str(wd(4)),'\pirad'];[num2str(angle(hd(4))),'\pirad']});
text(wd(5)-0.003,angle(hd(5))-1.3,...
{['\omega_p=',num2str(wd(5)),'\pirad'];[num2str(angle(hd(5))),'\pirad']});
hx3=xlabel('归一化数字角频率/\pirad');
hy2=ylabel('相值/\pirad');
ht3=title('相频特性');
ha3=gca;
gridon
%滤波
data=filter(b3,a3,Sig);
figure(3)
h4=plot(t,data,'k',t,Sig,'b',t,s_need,'r');
ht4=title('滤波效果');
xlabel('t/s','fontsize',s_xy);
ylabel('幅值','fontsize',s_xy);
legend('滤波后波形','滤波前波形','50Hz信号')
ha4=gca;
gridon
P_hou=fftlzr(data);
figure(4)
subplot(211)
h5=stem(2*pi/(N*Ts).*(0:
N-1),abs(P_qian));
hx4=xlabel('模拟角频率/rad/s');
ht5=title('滤波前幅度频谱');
axis([0,1.5e3,0,2e3])
text(f1*2*pi-20,abs(P_qian(f1*N*Ts+1))+200,[num2str(f1)'Hz']);
text(f2*2*pi-20,abs(P_qian(f2*N*Ts+1))+200,[num2str(f2)'Hz']);
text(f3*2*pi-20,abs(P_qian(f3*N*Ts+1))+200,[num2str(f3)'Hz']);
ha5=gca;gridon
subplot(212)
h6=stem(2*pi/(N*Ts).*(0:
N-1),abs(P_hou));
hx5=xlabel('模拟角频率/rad/s');
ht6=title('滤波后幅度频谱');
axis([0,1.5e3,0,2e3])
text(f1*2*pi-20,abs(P_hou(f1*N*Ts+1))+200,[num2str(f1)'Hz']);
text(f2*2*pi-20,abs(P_hou(f2*N*Ts+1))+200,[num2str(f2)'Hz']);
text(f3*2*pi-20,abs(P_hou(f3*N*Ts+1))+200,[num2str(f3)'Hz']);
ha6=gca;gridon
figure(5)
h7=plot(t,Sig,t,s_need,'r',t,s_1,'g',t,s_2,'c');
hx6=xlabel('t/s');
hy3=xlabel('幅值');
ht7=title('生成信号');
set(h7,'linewidth',2)
legend('合成信号','50Hz','100Hz','200Hz');
ha7=gca;gridon
%图像处理
h=[h1,h2,h3,h4',h5,h6,hsf,hsx];
ht=[ht1,ht2,ht3,ht4,ht5,ht6,ht7,htsf,htsx];
hx=[hx1,hx2,hx3,hx4,hx5,hx6,hxsf,hxsx];
hy=[hy1,hy2,hy3,hysf,hysx];ha=[ha1,ha2,ha3,ha4,ha5,ha6,ha7,hasf,hasx];
set(h([2:
6,9:
10]),'Linewidth',w_f)
set(ha,'Linewidth',w_axi,'Fontsize',s_f);
set([hxhyht],'Fontsize',s_xy,'FontWeight','bold');
set((1:
6),'outerposition',get(0,'screensize'),'color','w')