数字滤波器设计（数字信号处理上机实验室）
最编程
2024-04-14 17:46:42
...
close all
clc
clear
df=1;N=256;
Ts=1/N/df;L=N*Ts;t=[0:Ts:L]; %采样时刻
%信号
S=3+0.6*sin(2*pi*90*t+pi/2)+ sin(2*pi*50*t-pi/6)+0.7*sin(2*pi*80*t-pi/6);
%显示原始信号
plot(S);
title('原始信号');grid on;
figure;
Y = fft(S,N); %做FFT变换
Ayy = (abs(Y)); %取模
Ayy=Ayy/(N/2); %换算成实际的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)/Ts/N; %换算成实际的频率值
plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果
title('幅度-频率曲线图'); grid on;
figure;
Pyy=[1:N/2];
for i=1:N/2
Pyy(i)=phase(Y(i)); %计算相位
Pyy(i)=Pyy(i)*180/pi; %换算为角度
end
plot(F(1:N/2),Pyy(1:N/2)); %显示相位图
title('相位-频率曲线图'); grid on;
figure;
%%%%%%%%%%%%%%%%%%%%%%%数字滤波器设计6步
%(1)给定的数字滤波器的设计指标
wp=2*pi*50*Ts; %滤波器的通带截止频率
ws=2*pi*80*Ts; %滤波器的阻带截止频率
Rp=1;As=15; %滤波器的通阻带衰减指标
ripple=10^(-Rp/20); %滤波器的通带衰减对应的幅度值
Attn=10^(-As/20); %滤波器的阻带衰减对应的幅度值
%(2)转换为模拟滤波器的技术指标
Fs=2000;T=1/Fs;Omgp=wp*Fs;Omgs=ws*Fs;
%(3)确定模拟滤波器的最小阶数N和截止频率；
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s') %计算阶数n和截止频率
%(4) 根据最小阶数N计算模拟低通原型滤波器的系统传递函数
[z0,p0,k0]=buttap(n); %设计归一化的巴特沃思模拟滤波器原型
[ba1,aa1]=zp2tf(z0,p0,k0); %由零极点得到传递函数
%(5)利用模拟域频率变换法将模拟低通原型滤波器转换为实际模拟滤波器的系统传递函数；
[ba,aa]=lp2lp(ba1,aa1,Omgc); %变换为模拟低通滤波器
%(6)用脉冲响应不变法将模拟滤波器转换为数字滤波器
[bd,ad]=impinvar(ba,aa,Fs) %脉冲响应不变法
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%数字滤波器特性
[H,w]=freqz(bd,ad); %求数字系统的频率特性
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1);plot(w/pi,abs(H));
ylabel('|H|');title('幅度响应');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,wp,ws,2*ws]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2);plot(w/pi,angle(H)/pi);
ylabel('\phi');title('相位响应');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,wp,ws,2*ws]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,3);plot(w/pi,dbH);title('幅度响应(dB)');
ylabel('dB');xlabel('频率(\pi)');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,wp,ws,12*ws]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid
subplot(2,2,4);zplane(bd,ad);axis([-1.1,1.1,-1.1,1.1]);title('零极点图');
S1=filter(bd,ad,S);
figure;
plot(S1);
title('原始信号');grid on;
figure;
Y = fft(S1,N); %做FFT变换
Ayy = (abs(Y)); %取模
Ayy=Ayy/(N/2); %换算成实际的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)/Ts/N; %换算成实际的频率值
plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果
title('幅度-频率曲线图'); grid on;
figure;
Pyy=[1:N/2];
for i=1:N/2
Pyy(i)=phase(Y(i)); %计算相位
Pyy(i)=Pyy(i)*180/pi; %换算为角度
end
plot(F(1:N/2),Pyy(1:N/2)); %显示相位图
title('相位-频率曲线图'); grid on;
上一篇： 分析和模拟 RC 电路的频率响应，频率选择性网络特性测试
下一篇： 幅频特性和相频特性的物理意义