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

2024-04-08 15:56:22

文章目录

一、实验目的
二、实验内容及步骤
三、回答思考题

一、实验目的

掌握用窗函数法设计 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;
%调用xt*生具有加性噪声的信号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;
%调用xt*生具有加性噪声的信号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('滤除噪声后的信号波形');

窗函数法：