Digital Signal Process – 数字滤波器的设计(实验二:IIR 和FIR数字滤波器的设计和实现(2学时))

数字信号处理的另一个重要应用是数字滤波器。数字滤波器是一个运算过程,将输入数列按既定的要求转换成输出数列。在数字信号的处理时只需利用数字相加、乘以常数和延时等运算,就可以完全达到传输特性的要求。数字滤波器分为无限冲激响应(IIR)和有限冲激响应(FIR)两大类。本次实验的内容是数字滤波器设计,要求如下:

采集一段声音信号,长度 >15秒,取样频率 > 10kHz,分别设计一个IIR和FIR的带通滤波器,指标为 wp1=3kHz,wp2=4kHz,ws1=2.7kHz,ws2=4.3kHz,Ap=0.1dB,As=70dB。设计这个滤波器的传输函数 H(z),分析

  1. IIR和FIR滤波器的系数差异;
  2. 针对IIR滤波器,采用Butterworth型;
  3. 针对FIR滤波器,采用不同窗函数(矩形窗、Hanning、Hamming 和 Blackman 窗)时原始和滤波后信号的幅频特性。

通过本次实验,掌握以下知识:

  1. 复习和巩固 IIR、FIR数字滤波器的基本概念;
  2. 掌握 IIR、FIR 数字滤波器的设计方法;
  3. 掌握 IIR、FIR 数字滤波器的实现结构;
  4. Gibbs 效应发生的原因和影响;
  5. 不同类型的窗函数对滤波效果的影响,以及窗函数的选择。

下面附录一些和本实验相关的MATLAB函数


和音频有关的MATLAB函数

audioread

读取音频文件

Read the data back into MATLAB using audioread.

[y,Fs] = audioread(‘handel.wav’);

Note: y存储的是音频数据,通常是矩阵格式,Fs是音频采样率。

参考文档:https://cn.mathworks.com/help/matlab/ref/audioread.html?searchHighlight=audioread&s_tid=doc_srchtitle

audiowrite

写音频文件

Create a WAVE file from the example file handel.mat, and read the file back into MATLAB®.

Write a WAVE (.wav) file in the current folder.

load handel.mat filename = 'handel.wav';
audiowrite(filename,y,Fs);

Note: y存储的是音频数据,通常是矩阵格式,Fs是音频采样率。

参考文档:https://cn.mathworks.com/help/matlab/ref/audiowrite.html?searchHighlight=audiowrite&s_tid=doc_srchtitle

设计IIR滤波器相关的函数

设计butterworth带通滤波器的步骤:

  1. 由所要求的数字带通滤波器的边缘频率求出对应的数字频率 。
  2. 对数字频率做预畸变 。
  3. 将带通滤波器指标转化为相应归一化低通技术指标。
  4. 求出模拟低通滤波器的阶次,利用函数[N,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)   注意:Wp,Ws应该为步骤3中的。
  5. 设计模拟低通原型滤波器,其调用格式是  [z,p,k]=buttap(N)。 N是欲设计的低通原型滤波器的阶次,z,p,k分别是设计出的的极点、零点及增益。
  6. 求模拟低通原型滤波器传输函数的分子分母系数, [b,a]=zp2tf(z,p,k)。
  7. 求出带通滤波器传输函数的分子、分数系数。[B,A]=lp2bp(b,a,Wo, Wn)。
  8. 利用bilinear函数求出所要求的数字butterworth带通滤波器的传输函数H(z)。
  9. 求频率响应,利用Freqz函数。

buttord

Butterworth filter order and cutoff frequency

buttord calculates the minimum order of a digital or analog Butterworth filter required to meet a set of filter design specifications.

[n,Wn] = buttord(Wp,Ws,Rp,Rs) returns the lowest order, n, of the digital Butterworth filter with no more than Rp dB of passband ripple and at least Rs dB of attenuation in the stopband. The scalar (or vector) of corresponding cutoff frequencies, Wn, is also returned. Use the output arguments n and Wn in butter.

buttord 既可以计算butterworth数字滤波器的阶数和截止频率,又可以计算butterworth模拟滤波器的阶数和截止频率。

[n,Wn] = buttord(Wp,Ws,Rp,Rs)      % 计算数字butterworth滤波器的阶数和截止频率
[n,Wn] = buttord(Wp,Ws,Rp,Rs,'s')  % 计算模拟butterworth滤波器的阶数和截止频率

在设计数字滤波器,buttord输入参数中的Wp,Ws的取值范围是[0,1],其单位是π rad/s,1对应数字频率中的π。假如某个低通滤波器的通带边缘频率是0.3π,阻带边缘频率是0.35π,那么应另Wp = 0.3, Ws = 0.35,即将数字频率除以一个π。

在设计模拟滤波器时,buttord输入参数中的Wp和Ws的单位应该是rad/s,其他相同。

如何确定buttord的输入参数来设计低通、带通、高通、带阻数字滤波器,请查看参考文档中的两个表格。

参考文档:http://cn.mathworks.com/help/signal/ref/buttord.html?searchHighlight=buttord&s_tid=doc_srchtitle

buttap

Butterworth filter prototype

[z,p,k] = buttap(n) returns the poles and gain of an order n Butterworth analog lowpass filter prototype. The function returns the poles in the length n column vector p and the gain in scalar k. is an empty matrix because there are no zeros

buttap的输入为buttord计算出来的滤波器阶数n。

参考文档:https://cn.mathworks.com/help/signal/ref/buttap.html?searchHighlight=buttap&s_tid=doc_srchtitle

zp2tf

Convert zero-pole-gain filter parameters to transfer function form

zp2tf forms transfer function polynomials from the zeros, poles, and gains of a system in factored form.

zp2tf根据零点、极点计算传输函数。

调用方式:

[b,a] = zp2tf(z,p,k)

参考文档:https://cn.mathworks.com/help/signal/ref/zp2tf.html?searchHighlight=zp2tf&s_tid=doc_srchtitle

lp2bp

Transform lowpass analog filters to bandpass

lp2bp transforms analog lowpass filter prototypes with a cutoff angular frequency of 1 rad/s into bandpass filters with desired bandwidth and center frequency. The transformation is one step in the digital filter design process for the butter, cheby1, cheby2, and ellip functions.

lp2bp can perform the transformation on two different linear system representations: transfer function form and state-space form. In both cases, the input system must be an analog filter prototype.

lp2bp将模拟低通滤波器转化成数字带通滤波器。

调用方式:

[bt,at] = lp2bp(b,a,Wo,Bw)

For a filter with lower band edge w1 and upper band edge w2, use Wo = sqrt(w1*w2) and Bw = w2-w1.

参考文档:https://cn.mathworks.com/help/signal/ref/lp2bp.html?searchHighlight=lp2bp&s_tid=doc_srchtitle

bilinear

Bilinear transformation method for analog-to-digital filter conversion

The bilinear transformation is a mathematical mapping of variables. In digital filtering, it is a standard method of mapping the s or analog plane into the z or digital plane. It transforms analog filters, designed using classical filter design techniques, into their discrete equivalents.

调用方式:

[numd,dend] = bilinear(num,den,fs)

输入参数中的num和den分别是lp2bp的输出参数bt和at,fs是系统采样率

参考文档:https://cn.mathworks.com/help/signal/ref/bilinear.html?searchHighlight=bilinear&s_tid=doc_srchtitle

 

freqz

Frequency response of digital filter

freqz(hfilt) uses FVTool to plot the magnitude and unwrapped phase of the frequency response of the filter hfilt. If hfilt is a vector of filters, freqz plots the magnitude response and phase for each filter in the vector.

使用类似于[]=freqz()的形式,不会画图,使用类似freqz()的形式,会画出滤波器的幅度谱和相位谱。

调用方式可为:

Freqz(b, a)

其中,b和a是bilinear 计算出来的结果。

参考文档:https://cn.mathworks.com/help/dsp/ref/freqz.html?searchHighlight=freqz&s_tid=doc_srchtitle

filter

1-D digital filter

此 MATLAB 函数 使用分别由分子和分母系数 b 和 a 定义的有理传递函数筛选输入数据 x。

y = filter(b,a,x)

其中,b和a 是butter计算出来的结果,x为待滤波的数据,y是滤波之后的数据。

参考文档:https://cn.mathworks.com/help/matlab/ref/filter.html?searchHighlight=filter&s_tid=doc_srchtitle

设计FIR滤波器相关的函数

窗函数:

矩形窗函数:w = boxcar(N)

输入参数:窗的长度N

输出参数:窗序列向量w

 

汉明窗函数:w = hamming(N)

输入参数:窗的长度N

输出参数:窗序列向量w

 

汉宁窗函数:w = hanning(N)

输入参数:窗的长度N

输出参数:窗序列向量w

 

布莱克曼函数:w = blackman(N)

输入参数:窗的长度N

输出参数:窗序列向量w

fir1

Window-based FIR filter design

b = fir1(n,Wn) uses a Hamming window to design an nth-order lowpass, bandpass, or multiband FIR filter with linear phase. The filter type depends on the number of elements of Wn.

b = fir1(n,Wn,ftype) designs a lowpass, highpass, bandpass, bandstop, or multiband filter, depending on the value of ftype and the number of elements of Wn.

b = fir1(___,window) designs the filter using the vector specified in window and any of the arguments from previous syntaxes.

在调用形式中,n为滤波器的阶数,Wn为截止频率,其取值范围也是[0,1],转换方法同buttord中频率的转换方法,ftype为滤波器类型,取值为lowpass, highpass, bandpass, bandstop, multiband。滤波器的类型除了和ftype有关之外,还和Wn的维数有关,请注意两者要匹配。

输出为n阶滤波器系数向量b,这是一个截止频率为Wn线性相位FIR数字滤波器,Wn是从0到1之间的数,1对应奈式频率。对于带通滤波器或带阻滤波器,可将Wn定义为包含通频带边带频率的一个二元素向量,对于阻带滤波器,在“ftype”中附加“stop”字样。

向量window表示设计FIR滤波器所采用的的窗函数类型,以列向量形式表示。向量window的长度必须为n+1,若window缺省,则fir1使用汉明窗。

如果想设计hanning窗,则

b = fir1(n, Wn, hanning(n + 1))

参考文档:https://cn.mathworks.com/help/signal/ref/fir1.html?searchHighlight=fir1&s_tid=doc_srchtitle

freqz

Frequency response of digital filter

该函数的使用方法与设计IIR时差不多,但也稍有不同。可按照如下调用形式画出数字滤波器的幅度谱和相位谱

 

freqz(b)
freqz(dfilt.dffir(b))

其中b为fir1计算出来的结果

filter

1-D digital filter

Filter同设计IIR一样的形式调用。

y = filter(b,a,x)

不过不同的是,b是filter计算出来的结果,而a确实常数1,原因参看参考文档。同样,x为待滤波的数据,y是滤波之后的数据。故,在设计FIR数字滤波器的调用形式为

y = filter(b,1,x)

参考文档:https://cn.mathworks.com/help/signal/ref/filter.html

 

设计带通butterworth数字滤波器的例子

参考《数字信号处理(第二版)》 科学出版社 (门爱东 编著)P197 例4.13

设计代码如下:

%% FILE INTRODUCTION 
% file name:    exampleButterworth.m
% funtion:      design butterworth digital bandpass filter
% author:       chengkai
% last date:    2016/12/23
% create date:  2016/12/22
% release:      Second Release

clc,clear
close all;
%% 数字滤波器参数
fp1 = 300;
fp2 = 400;
fs1 = 200;
fs2 = 500;
butterAp = 3;       % dB
butterAs = 18;      % dB
fs = 2000;

%% 模拟频率转换成数字频率
wp1 = 2 * pi * fp1 / fs;
wp2 = 2 * pi * fp2 / fs;
ws1 = 2 * pi * fs1 / fs;
ws2 = 2 * pi * fs2 / fs;
%% IIR - Butterworth 滤波器的设计过程
omegaP=[2*fs*tan(wp1/2) 2*fs*tan(wp2/2)];
omegaS=[2*fs*tan(ws1/2) 2*fs*tan(ws2/2)];  
butterB = omegaP(2) - omegaP(1);
omega02 = omegaP(1) * omegaP(2);
if(omega02 ~= omegaS(1) * omegaS(2))
    % 调整其中一个阻带截止频率
    omegaS(1) = omega02 / omegaS(2);
end
% 将带通滤波器指标转化为相应归一化低通技术指标
lambdaP = 1;
lambdaS = (omegaS(2) - omegaS(1)) / (omegaP(2) - omegaP(1));
[N,omegaC]=buttord(lambdaP,lambdaS,butterAp,butterAs,'s');
[z,p,k]=buttap(N);  
[b,a]=zp2tf(z,p,k);  
[d,c]=lp2bp(b,a,sqrt(omega02), (butterB));
[bz,az]=bilinear(d,c,fs);  
freqz(bz,az);

代码和结果请点击进入:http://www.wydfx.xyz/wp-dsp/exampleButterworth.html