其他分享
首页 > 其他分享> > 基于MATLAB的数字信号处理(5) FIR数字滤波器设计及软件实现

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

作者：互联网

文章目录

一、实验目的

二、实验内容及步骤

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('幅度')

分析：

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

程序框图如下：

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]);      %设置刻度  右边

窗函数法：

等波纹最佳逼近法：

两种设计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/
本文仅用于交流学习，未经作者允许，禁止转载，更勿做其他用途，违者必究。
觉得文章对你有帮助、让你有所收获的话，期待你的点赞呀，不足之处，也可以在评论区多多指正。

标签：FIR,滤波器,Fs,数字,函数,阻带,设计,xt,MATLAB
来源： https://blog.csdn.net/fyfugoyfa/article/details/112924193