FIR滤波器和IIR滤波器

Alleyf

已于 2023-11-04 16:49:43 修改

阅读量5.2k

点赞数 9

文章标签： dsp开发

于 2022-10-14 00:09:00 首次发布

本文链接：https://blog.csdn.net/qq_61673686/article/details/127311983

版权

信号的时域处理

一、要求

选择子作业１中的音频信号，自行给定滤波器的系统函数，分别采用时域线性卷积和差分方程两种方法对音频信号进行滤波处理，比较滤波前后信号的波形和回放的效果。

二、实现思路

1.时域卷积法

分析原始音频信号的频谱，观察有用信号和噪声的分布，选择滤波器类型（低通、高通、带通、带阻等）
设定幅频响应下降三分贝的截止频率和滤波器阶数
选择合适的窗函数（包括Rectangular、Hanning、Hamming、Bartlett、Blackman、Kaiser、Gaussian、Flat-top等）
计算滤波器系数（滤波器的单位脉冲响应）
将滤波器系数与原始信号卷积进行滤波

2.差分方程法

分析原始音频信号的频谱，观察有用信号和噪声的分布，选择滤波器类型（低通、高通、带通、带阻等）
设置通带截止频率和阻带截止频率并归一化，设定通带纹波系数和阻带衰减系数
求出切比雪夫滤波器的阶数和下降3分贝通带截止频率
计算滤波器的差分方程系数
调用filter函数对原始信号进行滤波

三、实现过程

滤波前后对比图

FIR hamming窗低通滤波器前后对比图

IIR 切比雪夫低通滤波器前后对比图

1.时域卷积

首先导入原始音频信号并画出时域图和幅频响应曲线，分析其频谱分布，代码如下所示。

%% 数据导入及参数设置
[x,Fs]=audioread('D:\数据\Carmen_overture_noisy_8k_9.5k.wav');
Nx=length(x);
tx=(0:Nx-1)/Fs;%时域范围
df=Fs/length(tx); %计算谱线间隔
fx=tx*df;%频域范围
X=FFT(x,Nx,Fs,tx);%快速傅里叶变换
fx=(fx*Fs)/1e3;%只显示正频谱
X=abs(X);
% 绘出时域图和幅频响应
figure;
subplot(221)
plot(tx,x); 
title("原信号音频时域图");
xlabel("t/s"); 
subplot(222)
% plot(fx*Fs/1e3,abs(X),"blue")
plot(fx(1,1:ceil(Nx/2)),X(1:ceil(Nx/2),1),'b')
title("原信号音频相频响应");
xlabel("f/Khz");

观察其幅频响应发现有用信号的频率在8Khz以下，以上则为能量较大的高频噪声，因此我采用低通滤波器对其进行滤波。

FIR滤波器

设FIR滤波器的单位冲激响应h(n)为一个N点序列，0≤n≤N-1，则滤波器的系统函数为：
$H(\mathrm{z})=\sum_{k=0}^{N-1}h(\mathrm{k})*\mathrm{z}^{-k}$
直接型FIR滤波器的差分方程为：
$\mathrm{y(n)}=\sum_{m=0}^{N-1}\mathrm{h(m)x(n-m)}$
由差分方程可知只需要求出FIR滤波器的单位冲激响应h(n)，设置滤波阶数和3db通带截止频率fc，采用汉林窗函数处理，再使用fir1函数将其与待滤波信号卷积即可完成滤波，具体实现代码如下。

%% hamming窗低通滤波器滤除高频噪声
fc=7.95e3; %下降3分贝截止频率
h_n=fir1(1000,fc*2/Fs,"low");
y=conv(h_n,x);

Ny=length(y);%滤波后信号的长度
ty=(0:Ny-1)/Fs;%时域范围
Y=FFT(y,Ny,Fs,ty);%快速傅里叶变换求频谱
df=Fs/length(ty); %计算谱线间隔
fy=ty*df;%频域范围
fy=(fy*Fs)/1e3;%只显示正频谱
Y=abs(Y);
% 绘出滤波后的时域图和幅频响应
subplot(223)
plot(ty,y,"g");
title("滤波后信号音频时域图");
xlabel("t/s"); 
subplot(224)
stem(fx(1,1:ceil(Ny/2)),Y(1:ceil(Ny/2),1),'c','.')
title("滤波后信号音频相频响应");
xlabel("f/Khz"); 
%% 试听及保存
% sound(y,Fs);
audiowrite("D:\数据\FIRfilter.wav",y,Fs);
%% 快速傅里叶变换函数
function FourierTransform = FFT(signal,N,fs,t)
Y=fft(signal,N);
df=fs/length(t); %计算谱线间隔
f=t*df;%频域范围
Ys=abs(Y); %幅度响应
FourierTransform=Ys;
end

FIR滤波器的单位冲激响应如下图所示：

通过时域卷积得到滤波后的结果如下图所示：

由图可知，发现经过1000阶的3分贝截止频率为7.95e3Khz的低通滤波器处理后，8Khz及以上的高频噪声基本被完全滤除，成功从原始音频信号中提取到有用信号。

2.差分方程

首先也是先对原始音频信号做傅里叶变换分析其频谱，频谱结果同上，这里不在赘述。

IIR切比雪夫低通滤波

IIR(infinite impulse response)滤波器，无限冲激响应滤波器,IIR滤波器定义如下：
$\sum_{l=0}^{N}a(l)y(n-l)= \sum _{k=0}^{M}b(k)x(n-k)$
令 $a (0) = 1, a (l), l = 1, 2, ... N$ 取反后得：
$\sum _{k=0}^{M}b(k)x(n-k)+ \sum _{l=1}^{N}a(l)y(n-l)$
IIR直接II型滤波器的流程框图如下所示：

根据流程框图可知，需要求得差分方程系数向量a和b，这里我采用常见的切比雪夫滤波器设定通带频率、阻带频率以及纹波系数和衰减系数计算出系数a和b，再利用filter函数传入系数a、b和待滤波信号即可完成滤波，采用如下代码实现差分方程法实现IIR切比雪夫低通滤波器。

%% IIR-Chebyshev低通滤波器滤除高频噪声
Wp=7.6e3*2/Fs;%通带频率
Ws=8e3*2/Fs;%阻带频率
[n,Wc]=cheb1ord(Wp,Ws,2,40);%计算阶数和3db通带截止频率
[b,a]=cheby1(n,3,Wc);%切比雪夫滤波器返回差分方程系数
freqz(b,a);%绘出滤波器频谱图（幅频响应和相频响应）
y=filter(b,a,x);%滤除高频噪声
Y=FFT(y,N,Fs,t);
Y=fftshift(Y);
Y=abs(Y);
% 绘出滤波后的时域图和幅频响应
subplot(223)
plot(t,y,'b');
title("滤波后信号音频时域图");
xlabel("t/s"); 
subplot(224)
% plot(f*Fs/1e3,Y,'k');
plot(f(1,ceil(N/2):end),Y(ceil(N/2):end,1),'g')
title("滤波后信号音频幅频响应");
xlabel("f/Khz");