数字信号处理实验:数字滤波器的设计与应用 |
您所在的位置:网站首页 › 滤波器的设计实验报告 › 数字信号处理实验:数字滤波器的设计与应用 |
一. 实验目的 1.掌握模拟滤波器的设计方法,以及脉冲响应不变法和双线性变换法设计IIR数字滤波 器的方法,针对实际信号能设计相应的 IIR 数字滤波器,并按要求进行滤波。 2.掌握用窗函数法设计FIR数字滤波器的方法,并通过实验了解各种窗函数对滤波特性的影响;针对已知多频率正弦信号的频谱或实际信号,能设计相应的FIR滤波器按要求进行滤波。 二. 实验原理 1.基于模拟滤波器的IIR的数字滤波器设计(双线性变换法): a.确定抽样频率T。双线性变换法中参数T的选择和最终设计出的数字滤波器无关,因此可以取实际关系中的值,有时为了设计简单,常取T=2; b.按照 c.按照模拟滤波器的技术指标 d.将模拟滤波器 2.FIR数字滤波器的设计(窗函数设计法): a.依据理想的频率响应函数 b.依据阻带最小衰减要求,结合书中的表格,选择窗函数类型;依据过渡带宽的要求确定窗的长度N; c.加窗处理,得到设计结果为 MATLAB中提供了利用函数fir1实现窗函数设计FIR数字滤波器的通用应用格式,即 其中, 三. 实验内容 (1) 用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为:在通带内频率低于0.2π 时,最大衰减小于1dB;在阻带内[0.3π,π]频率区间上,最小衰减大于15dB。其中采样间隔为0.02π 。 (2) 用步骤(1)所设计的滤波器对实际心电图信号采样序列进行仿真滤波处理,观察总结滤波作用与效果,心电图信号为 xn={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-26,12,8,0,-16, -38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0, -2,-2,0,0,-2,-2,-2,-2,0} (3)已知信号 四. 实验步骤 ①用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为:在通带内频率低于0.2π 时,最大衰减小于1dB;在阻带内[0.3π,π] 频率区间上,最小衰减大于15dB。其中采样间隔为0.02π 。 ②用步骤①所设计的滤波器对实际心电图信号采样序列进行仿真滤波处理,观察总结滤波作用与效果,心电图信号为 xn={-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-26,12,8,0,-16, -38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0, -2,-2,0,0,-2,-2,-2,-2,0} 实验代码: clc,clear; %用双线性变换法设计一个巴特沃斯低通IIR数字滤波器 wp=0.2*pi;ws=0.3*pi; %数字滤波器截止频率 Ap=1;As=15; %衰减设置 T=0.02*pi;fs=1/T; %抽样间隔与抽样频率 Wp=(2/T)*tan(wp/2);Ws=(2/T)*tan(ws/2); %转换为模拟滤波器截止频率 [N,Wc]=buttord(Wp,Ws,Ap,As,'s'); %计算阶数和截止频率 [z,p,k]=buttap(N); %归一化原型滤波器设计 B=k*real(poly(z)); %分子多项式系数 A=real(poly(p)); %分母多项式系数 [Bs,As]=lp2lp(B,A,Wc); %去归一化得到模拟低通滤波器 [Bz,Az]=bilinear(Bs,As,fs); %数字低通滤波器系数 [Hz,w]=freqz(Bz,Az); %数字低通滤波器的频率响应 dbHz=20*log10(abs(Hz)/max(abs(Hz))); %化为分贝值 subplot(1,3,1);plot(w/pi,abs(Hz));grid on; set(gca,'xtick',[0 0.2 0.3 1]);set(gca,'xticklabel',[0 0.2 0.3 1]); set(gca,'ytick',[0 0.178 0.891 1]);set(gca,'yticklabel',[0 0.178 0.891 1]); xlabel('\omega/\pi');ylabel('|H(e^j^\omega)|'); subplot(1,3,2);plot(w/pi,angle(Hz));grid on; set(gca,'xtick',[0 0.2 0.3 1]);set(gca,'xticklabel',[0 0.2 0.3 1]); xlabel('\omega/\pi');ylabel('相位'); subplot(1,3,3);plot(w/pi,dbHz);grid on; axis([0,1,-80,5]); set(gca,'xtick',[0 0.2 0.3 1]);set(gca,'xticklabel',[0 0.2 0.3 1]); set(gca,'ytick',[-80 -15 -1 0]);set(gca,'yticklabel',[-80 -15 -1 0]); xlabel('\omega/\pi');ylabel('幅度(bB)'); % 对心电图数字信号进行滤波 x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,... -60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,... -4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0]; y=filter(Bz,Az,x); %滤波函数进行滤波 subplot(2,1,1);plot(x),axis([0 60 -100 40]),title('原始信号'); subplot(2,1,2);plot(y),axis([0 60 -100 40]),title('滤波后信号');仿真结果: 由图可知设计出的数字低通滤波器阻带满足指标要求,通带指标有富余,没有频谱混叠。但数字滤波器和模拟滤波器的幅频响应曲线形状有很大区别,这是频率非线性畸变引起的。
将离散的心电图转变为连续信号进行观察,从图中可以清晰的看出:原始信号的波形波动较大,存在或多或少的“突起”,但是经过低通滤波器后,过滤掉了阻带不满足要求的信号,滤波后的信号看起来更为的平滑,可见达到了滤除无效部分的效果。 ③已知信号: 低通滤波器实验代码: clc,clear; % 用汉明窗设计FIR低通滤波器 wc=0.35; %归一化6dB截止频率,相当于除以pi alpha=20; N=2*alpha+1; %FIR滤波器长度 hn=fir1(N-1,wc,'low',hamming(N)); %用汉明窗设计低通滤波器 omega=linspace(0,pi,512); %频率抽样512个点 mag=20*log10(abs(freqz(hn,1,omega))); %计算幅度频率响应的频率响应 subplot(321);stem((0:N-1),hn,'.');grid on; xlabel('n');ylabel('h(n)');title('单位抽样响应'); subplot(322);plot(omega/pi,mag);grid on; axis([0 1 -80 10]); set(gca,'xtick',[0 0.5 1]);set(gca,'xticklabel',[0 0.5 1]); set(gca,'ytick',[-80 -53 -6 0]);set(gca,'yticklabel',[-80 -53 -6 0]); xlabel('\omega/\pi');ylabel('dB');title('幅度频率响应'); % 采用该FIR低通滤波器对已知信号进行信号提取 a=42;b=1; fs=20000; %采样频率 f1=a*b*50;f2=a*b*100;f3=a*b*150; %待滤波余弦信号频率 t=(0:400)/fs; %定义时间步长 s=cos(2*f1*pi*t)+cos(2*f2*pi*t)+cos(2*f3*pi*t); sf=filter(hn,1,s); %使用filter函数对信号进行滤波 subplot(323);plot(t,s); %滤波前的信号图像 axis([0 0.005 -1.5 3.5]); xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图'); Fs=fft(s,512);AFs=abs(Fs);f=fs/512*(0:255); subplot(324);plot(f,AFs(1:256)); %滤波前的信号频谱图 xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图'); subplot(325);plot(t,sf); %滤波后的信号图像 axis([0 0.005 -1.1 1.1]); xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图'); Fsf=fft(sf,512); %滤波后的信号频谱图 AFsf=abs(Fsf); %信号频谱图的幅值 subplot(326);plot(f,AFsf(1:256)); %滤波后的信号频谱图 xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');仿真结果: 带通滤波器实验代码: clc,clear; % 用汉明窗设计FIR带通滤波器 wc1=0.35;wc2=0.5; %归一化6dB截止频率,相当于除以pi alpha=20; N=2*alpha+1; %FIR滤波器长度 hn=fir1(N-1,[wc1 wc2],'bandpass',hamming(N)); %用汉明窗设计带通滤波器 omega=linspace(0,pi,512); %频率抽样512个点 mag=20*log10(abs(freqz(hn,1,omega))); %计算幅度频率响应的频率响应 subplot(321);stem((0:N-1),hn,'.');grid on; xlabel('n');ylabel('h(n)');title('单位抽样响应'); subplot(322);plot(omega/pi,mag);grid on; axis([0 1 -80 10]); set(gca,'xtick',[0 0.5 1]);set(gca,'xticklabel',[0 0.5 1]); set(gca,'ytick',[-80 -53 -6 0]);set(gca,'yticklabel',[-80 -53 -6 0]); xlabel('\omega/\pi');ylabel('dB');title('幅度频率响应'); % 采用该FIR带通滤波器对已知信号进行信号提取 a=42;b=1; fs=20000; %采样频率 f1=a*b*50;f2=a*b*100;f3=a*b*150; %待滤波余弦信号频率 t=(0:400)/fs; %定义时间步长 s=cos(2*f1*pi*t)+cos(2*f2*pi*t)+cos(2*f3*pi*t); sf=filter(hn,1,s); %使用filter函数对信号进行滤波 subplot(323);plot(t,s); %滤波前的信号图像 axis([0 0.005 -1.5 3.5]); xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图'); Fs=fft(s,512);AFs=abs(Fs);f=fs/512*(0:255); subplot(324);plot(f,AFs(1:256)); %滤波前的信号频谱图 xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图'); subplot(325);plot(t,sf); %滤波后的信号图像 axis([0 0.005 -1.1 1.1]); xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图'); Fsf=fft(sf,512); %滤波后的信号频谱图 AFsf=abs(Fsf); %信号频谱图的幅值 subplot(326);plot(f,AFsf(1:256)); %滤波后的信号频谱图 xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');仿真结果: 高通滤波器实验代码: clc,clear; % 用汉明窗设计FIR高通滤波器 wc=0.5; %归一化6dB截止频率,相当于除以pi alpha=20; N=2*alpha+1; %FIR滤波器长度 hn=fir1(N-1,wc,'high',hamming(N)); %用汉明窗设计高通滤波器 omega=linspace(0,pi,512); %频率抽样512个点 mag=20*log10(abs(freqz(hn,1,omega))); %计算幅度频率响应的频率响应 subplot(321);stem((0:N-1),hn,'.');grid on; xlabel('n');ylabel('h(n)');title('单位抽样响应'); subplot(322);plot(omega/pi,mag);grid on; axis([0 1 -80 10]); set(gca,'xtick',[0 0.5 1]);set(gca,'xticklabel',[0 0.5 1]); set(gca,'ytick',[-80 -53 -6 0]);set(gca,'yticklabel',[-80 -53 -6 0]); xlabel('\omega/\pi');ylabel('dB');title('幅度频率响应'); % 采用该FIR高通滤波器对已知信号进行信号提取 a=42;b=1; fs=20000; %采样频率 f1=a*b*50;f2=a*b*100;f3=a*b*150; %待滤波余弦信号频率 t=(0:400)/fs; %定义时间步长 s=cos(2*f1*pi*t)+cos(2*f2*pi*t)+cos(2*f3*pi*t); sf=filter(hn,1,s); %使用filter函数对信号进行滤波 subplot(323);plot(t,s); %滤波前的信号图像 axis([0 0.005 -1.5 3.5]); xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图'); Fs=fft(s,512);AFs=abs(Fs);f=fs/512*(0:255); subplot(324);plot(f,AFs(1:256)); %滤波前的信号频谱图 xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图'); subplot(325);plot(t,sf); %滤波后的信号图像 axis([0 0.005 -1.1 1.5]); xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图'); Fsf=fft(sf,512); %滤波后的信号频谱图 AFsf=abs(Fsf); %信号频谱图的幅值 subplot(326);plot(f,AFsf(1:256)); %滤波后的信号频谱图 xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');仿真结果: 带阻滤波器实验代码: clc,clear; % 用汉明窗设计FIR带阻滤波器 wc1=0.3;wc2=0.55; %归一化6dB截止频率,相当于除以pi alpha=20; N=2*alpha+1; %FIR滤波器长度 hn=fir1(N-1,[wc1 wc2],'stop',hamming(N)); %用汉明窗设计带阻滤波器 omega=linspace(0,pi,512); %频率抽样512个点 mag=20*log10(abs(freqz(hn,1,omega))); %计算幅度频率响应的频率响应 subplot(321);stem((0:N-1),hn,'.');grid on; xlabel('n');ylabel('h(n)');title('单位抽样响应'); subplot(322);plot(omega/pi,mag);grid on; axis([0 1 -80 10]); set(gca,'xtick',[0 0.5 1]);set(gca,'xticklabel',[0 0.5 1]); set(gca,'ytick',[-80 -53 -6 0]);set(gca,'yticklabel',[-80 -53 -6 0]); xlabel('\omega/\pi');ylabel('dB');title('幅度频率响应'); % 采用该FIR带阻滤波器对已知信号进行信号提取 a=42;b=1; fs=20000; %采样频率 f1=a*b*50;f2=a*b*100;f3=a*b*150; %待滤波余弦信号频率 t=(0:400)/fs; %定义时间步长 s=cos(2*f1*pi*t)+cos(2*f2*pi*t)+cos(2*f3*pi*t); sf=filter(hn,1,s); %使用filter函数对信号进行滤波 subplot(323);plot(t,s); %滤波前的信号图像 axis([0 0.005 -1.5 3.5]); xlabel('时间/秒');ylabel('幅度');title('信号滤波前时域图'); Fs=fft(s,512);AFs=abs(Fs);f=fs/512*(0:255); subplot(324);plot(f,AFs(1:256)); %滤波前的信号频谱图 xlabel('频率/赫兹');ylabel('幅度');title('信号滤波前频域图'); subplot(325);plot(t,sf); %滤波后的信号图像 axis([0 0.005 -2 2.2]); xlabel('时间/秒');ylabel('幅度');title('信号滤波后时域图'); Fsf=fft(sf,512); %滤波后的信号频谱图 AFsf=abs(Fsf); %信号频谱图的幅值 subplot(326);plot(f,AFsf(1:256)); %滤波后的信号频谱图 xlabel('频率/赫兹');ylabel('幅度');title('信号滤波后频域图');仿真结果: 总结:使用窗函数是为了消除无限序列的截短而引起的吉布效应。加窗后,使频响产生一过渡带,其宽度正好等于窗的频响的主瓣宽度。出现肩峰时,肩峰两侧形成起伏振荡,其振荡幅度取决于旁瓣的相对幅度。 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |