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

叶庭云

发布于 2021-12-01 07:12:56

9.5K00

代码可运行

文章被收录于专栏：Python进阶之路Python进阶之路

运行总次数：0

代码可运行

文章目录

一、实验目的

掌握用窗函数法设计 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('滤除噪声后的信号波形');