基于MATLAB的数字信号处理(5) FIR数字滤波器设计及软件实现

一、实验目的

  • 掌握用窗函数法设计 FIR 数字滤波器的原理和方法
  • 掌握用等波纹最佳逼近法设计 FIR 数字滤波器的原理和方法
  • 掌握 FIR 滤波器的快速卷积实现原理
  • 学会调用 MATLAB 函数设计与实现 FIR 滤波器

二、实验内容及步骤

  • 认真复习用窗函数法和等波纹最佳逼近法设计 FIR 数字滤波器的原理;
  • 调用信号产生函数 xtg 产生具有加性噪声的信号x(t),显示 x(t) 及其频谱,如下图所示:
function xt=xtg(N)

%信号x(t)产生,并显示信号的幅频特性曲线
%xt=xtg(N) 产生一个长度为N,有加性高频噪声的单频调幅信号x(t),采样频率Fs=1000Hz
%载波频率fc=Fs/10=100Hz,调制正弦波频率f0=fc/10=10Hz.
N=1000;
Fs=1000;T=1/Fs;Tp=N*T;
t=0:T:(N-1)*T;
fc=Fs/10;f0=fc/10;  %载波频率fc=Fs/10,单频调制信号频率为f0=Fc/10;
mt=cos(2*pi*f0*t);  %产生单频正弦波调制信号mt,频率为f0
ct=cos(2*pi*fc*t);  %产生载波正弦波信号ct,频率为fc
xt=mt.*ct;          %相乘产生单频调制信号xt
nt=2*rand(1,N)-1;   %产生随机噪声nt

%======= 设计高通滤波器hn,用于滤除噪声nt中的低频成分,生成高通噪声 =======
fp=150; fs=200;Rp=0.1;As=70;  %滤波器指标
fb=[fp,fs];m=[0,1]; 		  %计算remezord函数所需参数f,m,dev
dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];
[n,fo,mo,W]=remezord(fb,m,dev,Fs);	%确定remez函数所需参数
hn=remez(n,fo,mo,W);	    %调用remez函数进行设计,用于滤除噪声nt中的低频成分
yt=filter(hn,1,10*nt);      %滤除随机噪声中低频成分,生成高通噪声yt

%噪声加信号
xt=xt+yt;           
fst=fft(xt,N);k=0:N-1;f=k/Tp;
subplot(211);plot(t,xt,'k');grid;xlabel('t/s');ylabel('x(t)');
axis([0,Tp/2,min(xt),max(xt)]);title('(a) 信号加噪声波形')
subplot(212);plot(f,abs(fst)/max(abs(fst)),'r');grid;title('(b) 信号加噪声频谱')
axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度')

  • 设计低通滤波器,从高频噪声中提取x(t)中的单频调幅信号,要求信号幅频失真小于0.1 dB,将噪声频谱衰减60 dB。先观察 x(t) 的频谱,确定滤波器指标参数。
  • 根据滤波器指标选择合适的窗函数,计算窗函数的长度N, 调用 MATLAB 函数 fir1 设计一个 FIR 低通滤波器,并编写程序,调用 MATLAB快速卷积函数 fftfilt 实现对 x(t) 的滤波。绘图显示滤波器的频率响应特性曲线、滤波器输出信号的幅频特性图和时域波形图。
  • 滤波器指标不变,但改用等波纹最佳逼近法,调用 MATLAB 函数 remezord 和 remez 设计 FIR 数字滤波器。并比较两种设计方法设计的滤波器阶数。

分析:

MATLAB函数 fir1 和 fftfilt 的功能及其调用格式可以用 help 命令查阅;

  • 取采样频率Fs =1000Hz,采样周期T=1/Fs ;
  • 根据实验要求,可选择滤波器指标参数:通带截止频率fp=120Hz,阻带截止频率fs=150Hz,换算成数字频率,通带截止频率ωp =2πfpT=0.24π rad,通带最大衰减为0.1dB,阻带截止频率ωs=2πfsT=0.3π rad,阻带最小衰减为60dB;

程序框图如下:

clear;
%调用xtg产生具有加性噪声的信号x(t)
%x(t)长度N=1000,并显示x(t)信号波形及其频谱
N=1000;xt=xtg(N);Fs=1000;T=1/Fs;Tp=N*T;
k=0:N-1;f1=k/Tp; 
 
%初始化指标
%通带截止频率fp 阻带截止频率fs 通带最大衰减 阻带最小衰减 采样频率Fs=1000Hz
fp=120;fs=150;Rp=0.2;As=60;Fs=1000;  
 
%用窗函数法设计滤波器
wc=(fp+fs)/Fs;      %理想低通滤波器截止频率(关于pi归一化)
B=2*pi*(fs-fp)/Fs;  %过渡带宽度指标
Nw=ceil(11*pi/B);   %blackman窗的长度N
hn=fir1(Nw-1,wc,blackman(Nw)); %阻带最小衰减为60dB 可以用blackman窗
disp(Nw-1);             %窗函数法设计滤波器 滤波器长度-1=阶数
Hw=abs(fft(hn,1024));   %求设计的滤波器频率特性
ywt=fftfilt(hn,xt,N);   %调用函数fftfilt对xt滤波
 
%窗函数法设计法的绘图部分(滤波器损耗函数, 滤波器输出信号波形)
f=(0:1023)*Fs/1024;
figure(2);subplot(3,1,1);
plot(f,20*log10(Hw/max(Hw)),'g');grid;title('窗函数法低通滤波器幅频特性')
axis([0,Fs/2,-120,20]);
xlabel('f/Hz');ylabel('幅度')
t=(0:N-1)/Fs;Tp=N/Fs; 
subplot(3,1,2);
plot(t,ywt,'r');grid;
axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_w(t)');
title('滤除噪声后的信号波形');
 
subplot(3,1,3);
fxt=fftshift(fft(ywt,1000));
stem(f1-500,abs(fxt)/max(abs(fxt)),'k.');grid;title('滤除噪声后的信号频谱');
axis([0,500,0,1]);      %设置刻度  右边
 
%等波纹最佳逼近法设计滤波器
fb=[fp,fs]; m=[1,0];        
%确定remezord函数所需参数f,m,dev
dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; %十进制表示的
[Ne,fo,mo,W]=remezord(fb,m,dev,Fs);
disp(Ne-1);   %波纹最佳逼近法设计滤波器 滤波器长度
%确定remez函数所需参数
hn=remez(Ne-1,fo,mo,W);  %调用remez函数进行设计
Hw=abs(fft(hn,1024));  %求设计的滤波器频率特性
yet=fftfilt(hn,xt,N);  %调用函数fftfilt对xt滤波
 
%等波纹设计法的绘图部分(滤波器损耗函数, 滤波器输出信号yw(t)波形)
f=(0:1023)*Fs/1024;
figure(3);subplot(3,1,1);
plot(f,20*log10(Hw/max(Hw)),'g');grid;title('等波纹法低通滤波器幅频特性');
axis([0,Fs/2,-80,10]);
xlabel('f/Hz');ylabel('幅度')
t=(0:N-1)/Fs;Tp=N/Fs;
subplot(3,1,2)
plot(t,yet,'r');grid;
axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_e(t)');
title('滤除噪声后的信号波形');
 
subplot(3,1,3);
fet=fftshift(fft(yet,1000));
stem(f1-500,abs(fet)/max(abs(fet)),'k.');grid;title('滤除噪声后的信号频谱');
axis([0,500,0,1]);      %设置刻度  右边

窗函数法:

等波纹最佳逼近法:

  • 选择滤波器指标参数:通带截止频率fp=120Hz,阻带截止频率fs=150 Hz。代入采样频率Fs=1000Hz,将模拟频率换算成数字频率,通带截止频率ωp=2πfpT=0.24π rad,通带最大衰减为0.1dB,阻带截止频ωs=2πfsT=0.3π rad,阻带最小衰减为60dB。blackman窗阻带最小衰减为74dB,而指标仅为60dB,所以可以选取 blackman 窗函数,与信号产生函数 xtg 相同,采样频率取Fs=1000Hz。
  • 用窗函数法设计滤波器,滤波器长度Nw=184,等波纹最佳逼近法设计滤波器,滤波器长度Ne=83。两种方法设计的滤波器都能有效地从噪声中提取信号,但等波纹最佳逼近法设计的滤波器阶数低得多。当然,滤波实现的运算量以及时延也小得多。

两种设计FIR滤波器的方法(窗函数法和等波纹最佳逼近法)的分析比较:

  • 用窗函数法设计的滤波器,如果在阻带截止频率附近刚好满足,则离开阻带截止频率越远,阻带衰减富裕量越大,即存在资源浪费;
  • 用等波纹最佳逼近法设计的滤波器,其通带和阻带均为等波纹特性,且通带最大衰减和阻带最小衰减可以分别控制,所以其指标均匀分布,没有资源浪费,所以其阶数低得多。而在窗函数设计法中,通带、阻带波纹的改善是以加宽过渡带为代价的。

窗函数设计法中还可以改用 Kaiser 窗,设置β参数为5.658,此时阻带最小衰减为60dB,刚好符合指标,MATLAB绘图如下:

clear;
%调用xtg产生具有加性噪声的信号x(t)
%x(t)长度N=1000,并显示x(t)信号波形及其频谱
N=1000;xt=xtg(N);Fs=1000;T=1/Fs;Tp=N*T;
k=0:N-1;f1=k/Tp;
%初始化指标
%通带截止频率fp 阻带截止频率fs 通带最大衰减 阻带最小衰减 采样频率Fs=1000Hz
fp=120;fs=150;Rp=0.2;As=60;Fs=1000;  

%用窗函数法设计滤波器
wc=(fp+fs)/Fs;        %理想低通滤波器截止频率(关于pi归一化)
B=2*pi*(fs-fp)/Fs;    %过渡带宽度指标
Nw=ceil(7.24*pi/B);   %kaiser窗的长度N
hn=fir1(Nw-1,wc,kaiser(Nw,5.658)); %阻带最小衰减为60dB 用kaiser窗 设置β参数
disp(Nw-1);             %窗函数法设计滤波器 滤波器长度-1=阶数
Hw=abs(fft(hn,1024));   %求设计的滤波器频率特性
ywt=fftfilt(hn,xt,N);   %调用函数fftfilt对xt滤波

%窗函数法设计法的绘图部分(滤波器损耗函数, 滤波器输出信号波形)
f=(0:1023)*Fs/1024;
figure(2);subplot(3,1,1);
plot(f,20*log10(Hw/max(Hw)),'g');grid;title('窗函数法低通滤波器幅频特性')
axis([0,Fs/2,-120,20]);
xlabel('f/Hz');ylabel('幅度')
t=(0:N-1)/Fs;Tp=N/Fs; 
subplot(3,1,2);
plot(t,ywt,'r');grid;
axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_w(t)');
title('滤除噪声后的信号波形');

subplot(3,1,3);
%fxt=fftshift(fft(hn,1000));
fxt=fftshift(fft(ywt,1000));
stem(f1-500,abs(fxt)/max(abs(fxt)),'k.');grid;title('滤除噪声后的信号频谱');
axis([0,500,0,1]);      %设置刻度  右边

%等波纹最佳逼近法设计滤波器
fb=[fp,fs]; m=[1,0];  		
%确定remezord函数所需参数f,m,dev
dev=[(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; %十进制表示的通带和阻带的波动
[Ne,fo,mo,W]=remezord(fb,m,dev,Fs);
disp(Ne-1);   %波纹最佳逼近法设计滤波器 滤波器长度

%确定remez函数所需参数
hn=remez(Ne,fo,mo,W);  %调用remez函数进行设计
Hw=abs(fft(hn,1024));  %求设计的滤波器频率特性
yet=fftfilt(hn,xt,N);  %调用函数fftfilt对xt滤波

%等波纹设计法的绘图部分(滤波器损耗函数, 滤波器输出信号yw(t)波形)
f=(0:1023)*Fs/1024;
figure(3);subplot(2,1,1);
plot(f,20*log10(Hw/max(Hw)),'g');grid;title('等波纹法低通滤波器幅频特性');
axis([0,Fs/2,-80,10]);
xlabel('f/Hz');ylabel('幅度')
t=(0:N-1)/Fs;Tp=N/Fs;
subplot(2,1,2);
plot(t,yet,'r');grid;
axis([0,Tp/2,-1,1]);xlabel('t/s');ylabel('y_e(t)');
title('滤除噪声后的信号波形');

窗函数法:

等波纹最佳逼近法:


此时滤波器长度Nw=121,相比 blackman 窗的滤波器长度,减少了很多。滤波实现的运算量以及时延也小了些。

三、回答思考题

如果给定通带截止频率和阻带截止频率以及阻带最小衰减,如何用窗函数法设计线性相位低通滤波器?请写出设计步骤:

解释为什么对同样的技术指标,用等波纹最佳逼近法设计的滤波器阶数低?

作者:叶庭云
CSDN:https://yetingyun.blog.csdn.net/
本文仅用于交流学习,未经作者允许,禁止转载,更勿做其他用途,违者必究。
觉得文章对你有帮助、让你有所收获的话,期待你的点赞呀,不足之处,也可以在评论区多多指正。

已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 书香水墨 设计师:CSDN官方博客 返回首页