【医学信号处理与MATLAB（2）】数字滤波器的设计、FIR/IIR滤波器

本节主要介绍的是数字滤波器的设计，以下是对应信号处理流程中的框图

目录

数字滤波器的设计(Digital filter design)

滤波的目的

• 将信号的频谱成分作调整，去除某些频率。

如下图中的EEG信号，左边是时域图像，右边是频域图像，EEG中我们通常感兴趣的频率范围在50Hz左右，所以我们需要一种滤波器，在我们不感兴趣的频带将信号滤除掉。
图2中红线为一个低通滤波器（可以通过滤波器与信号点乘实现滤波）。而图3是低通滤波之后的EEG信号，图4为滤波后信号经过FFT的结果，可以看出基本只保留了0-50Hz的部分

• 移除高频仪器或环境噪音
• 移除低频信号漂移成分

我们通常会利用低通滤波器来去除这部分的噪音。如下图中的例子：

上图中的信号存在信号漂移等现象，这通常是很讨厌的，因为我们在事件侦测的时候不希望这种现象出现，所以我们希望消除这种现象，便使用高通滤波器（如图3，一般设置为1Hz），将这二者结合，便可以使用带通滤波实现

傅里叶变换的特性

傅里叶变换就是将信号利用cos()和sin()作为基底进行一个展开的动作，无论是周期性还是非周期信号，都可以使用。公式如下所示：

傅里叶转换有一些非常好用的特性，如下：

卷积(Convolution)：直观上来说就是不停做一个叠加的动作，类似于for-loop，而时域上的卷积，到频域上就变成了相乘的动作，这与滤波的原理是非常类似的
平移特性(Time shifting)：这个就比较简单了

Z变换的特性

如果说傅里叶变换是把信号投影到cos和sin上，那么Z变换就是把信号投影到z域上，其公式如下：

Z变换的特性与傅里叶变换相同

滤波器的运算原理

滤波器的主要运算原理如下图所示：

其中x(n)是离散原始信号，n=1:N，在时域上，加滤波器之后我们使用的是滤波器与信号的卷积，这通常是比较难以解决的，所以我们可以借助Z变换的特性，将时域原始信号转换成Z域原始信号，这样的话我们只需要计算滤波器与原始信号的点乘。
所以这也告诉我们，可以在Z域条件下设计滤波器，然后再还原到时域上，这里面便用到了我们的老朋友——离散传递函数。

离散传递函数

离散传函的定义如下图所示：
（这一块忘了的话去翻翻胡寿松的自控吧= =貌似笔者也忘了）

其中H(z)称为频率响应（单位冲激函数的情况下），而Y(z)就是我们滤波后的信号。

输入输出函数

理论上，我们使用al和bk就能决定一个滤波器。输出时域信号包含两种成分：原始信号和输出前面的信号(e.g. y(9) y(8)…)
滤波器的阶数为K和L，代表了al和bk的个数，也决定了滤波器的长度。理论来说，数字越大，会趋近与理想状态下的滤波器。
关于上述式子证明如下：(有兴趣可以看看)

传递函数(滤波器)的频谱分析

我们可以先看傅里叶变换的传函作为一个例子：

在频域上，我们如果想看H(f)的频率成分，我们可以经过代换后，分别对Y(f)和X(f)中的系数序列bk和al做傅里叶变换，再把它们点除之后得到最终的滤波器。
再给一个例子，请画出下列传递函数的频谱和相位变化：

我们可以在上面分别找出系数序列，然后分别进行傅里叶变换，代码如下所示：

clear,close all

%% initialize parameters
samplerate=1000;
N=512; %采样点数，必须是偶数，最好为2的幂次方

%% 定义传递函数H中的系数a和b
a=[1 -0.2 0.8];
b=[0.2 0.5];

%% 选择1：利用fft计算频谱H
H=fft(a,N)./fft(b,N);  %注意这里用点除，元素相除
mag=20log10(abs(H));  %取abs是因为H的值是一个复数，有时部和虚部，幅度的dB值
phase=angle(H)*2*pi;  %相位,直接用MATLAB中的angle函数，主要是看相位漂移
faxis=samplerate/2*linspace(0,1,N/2) %频率轴，注意写法

%% 绘制出频谱图
figure,
subplot(2,1,1),plot(faxis,mag(1:N/2))
xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
grid on

subplot(2,1,2),plot(faxis,phase(1:N/2))
xlabel('Frequency (Hz)'),ylabel('Phase (deg.)')
grid on

执行之后，绘制出来的曲线如下图所示：

这个幅频特性曲线有点像我们认知中的带通滤波器，但是不是理想的，因为阶数非常低，衰减情况可以假设图上一个点为-10dB，那么衰减倍数=10^(-10/20)=0.3162倍。

另一种方式：（直接用matlab函数实现）

[H,faxis]=freqz(b,z,N,'whole',samplerate);

它的整个代码如下所示：

clear,close all

%% initialize parameters
samplerate=1000;
N=512; %采样点数，必须是偶数，最好为2的幂次方

%% 定义传递函数H中的系数a和b
a=[1 -0.2 0.8];
b=[0.2 0.5];

%% 选择2，使用freqz函数计算频谱
[H,faxis]=freqz(b,z,N,'whole',samplerate);
mag=20log10(abs(H));
phase=angle(H)*2*pi;
faxis=samplerate/2*linear(0,1,N/2);
figure,
subplot(2,1,1),plot(faxis,mag(1:N/2));
xlabel('Frequency (Hz)'),ylabel('Magnitude (dB)')
grid on

subplot(2,1,2),plot(faxis,phase(1:N/2),r);
xlabel('Frequency (Hz)'),ylabel('Phase (deg.)')
grid on

有限脉冲响应(FIR)滤波器

FIR的传递函数和输入输出

FIR=finite impluse response
有了上面的铺垫，就可以很自然地引出FIR了
下面是FIR的传递函数与输入输出函数：

与上面对比一下可以发现，它其实是把a系数的影响给设为了1，只考虑了b序列的影响。
FIR的优点有：

• 稳定（因为分母不会为0）
• 线性相位变化

缺点：

• 需要较高阶数来满足滤波条件
• 计算效率较低

FIR滤波器设计

1. 将设定好的频率响应，反推算出传递函数的系数b
2. 系数的计算，可以通过对想得到的频率响应做反傅里叶变换取得（频域->时域），如下图所示：（前提是时域上sinc函数的长度为无限长，结果才是一个方波）
   

其中，sinc函数的公式为：

所以，我们可以利用sinc函数设计几种我们经常用到的FIR滤波器

fc是频域响应中我们想要切断的频率数（分为最大和最小）
Ts是采样周期
L是滤波器长度

如果只使用FIR的话，可能导致的结果为：

可以看出，上述幅频曲线里面有一些ripple，而这也是FIR的一个特性，我们可以使用窗函数来解决这样的问题，如下：

比较下面两者的差异：

H=fft(b,N);
HW=fft(b.*hamming(L+1）',N)

关于加上窗函数的FIR，MATLAB程序如下：

%% initial parameters
samplerate=400;
N=512;   %512个点
L=128;   %滤波器的长度
fl=60;   %low截止频率
fh=120;   %high截止频率
Ts=1/samplerate;    %采样周期

%%带通滤波器（其它三种同理，可以自己写写看）
b=[];
for k=1:L+1   %这里加1无所谓
   n=k-L/2;   %滤波器括号里面的数值，公式参考上面
   if n==0    %分母为零，就不能算了
      b(k)=fh-fl;   
   else
      b(k)=sin(2*pi*fh*Ts*n)/(pi*n)-sin(2*pi*fl*Ts*n)/(pi*n)

%% 计算幅频特性曲线with/without窗函数
H=fft(b,N);   %without窗函数
mag=20*log10(abs(H));

HW=fft(b.*hamming(L+1)',N);    %with窗函数
magW=20*log10(abs(HW));

faxis=samplerate/2*linspace(0,1,N/2);  %频率轴

MATLAB滤波器设计函数

对于ppt中下面那句话的理解：假如滤波器的长度为128，那么数据的长度至少为128*3(filtfilt函数)

MATLAB中定义函数的规则

function y = mean(x,dim)

！！要注意，这个函数的名称一定要与函数所在的.m文件名称一致，这非常重要，因为MATLAB在调用函数的时候是看文件名的。
可以在function下面、代码上面写一些注释，这样在命令行输入help的时候可以弹出来这些注释。

FIR滤波器函数实例

1、3阶段FIR
现在已经有了写好的FIR/IIR滤波器，附件里有，举个例子

function [sigfilter,n]=filter_3sFIR(sig,f,a,dev,fs)

sig-原始信号
f-截止频率
a-这个band是要还是不要，要的话1，不要的话0 dev-容忍ripple的程度
fs-采样频率 n-阶数

后面把ppt放上来更好理解了

注意下图中f,dev和a三个参数怎么设定

2、2阶段FIR

sigfilter=filter_2sFIR(sig,f,fs,n,type,win)

fs-采样率samplerate
n-滤波器的阶数（FIR可能要100+阶）
type-lowpass还是highpass(‘low’or’high’)
win-是否加窗（例如汉明窗）

无限脉冲响应(IIR)滤波器

IIR的传递函数和输入输出

IIR的传递函数如下图所示：

比FIR多了分母al系数
输入输出函数：（ppt截的，有点乱）

IIR的优点有：

• 仅需较低阶数来满足滤波条件
• 计算效率较高，但运算要持续用到整段信号

缺点：

• 可能不稳定，分母项为0
• 非线性相位变化（相频曲线中间不是斜线，可能是个曲线）

同阶数FIR与IIR比较

理想状态下应该是图中的方波，而明显能够看出，在阶数相同的时候，IIR的效果要比FIR好

MATLAB滤波器设计函数

巴特沃斯、切比雪夫、椭圆滤波器balabala

IIR滤波器函数实例

%% initialize parameters
samplerate=1000; % in Hz
N=1024; % data length

sinefreq1=50; % in Hz
sinefreq2=100; % in Hz
sinefreq3=200; % in Hz

fl=75;   % low-cutoff frequency
fh=165;  % high-cutoff frequency
trans_width=20; % in Hz. It is a half of transition band. if data length is not long enough, increase trans_width

rp=1;  % in dB
rs=40; % in dB

%% generate simulated signals 
t=[1:N]/samplerate;
sig1=sin(2*pi*sinefreq1*t);
sig2=sin(2*pi*sinefreq2*t);
sig3=sin(2*pi*sinefreq3*t);

data=sig1+sig2+sig3;

%% apply 3-stage IIR filter processing
[data_3sIIR,forder] = filter_3sIIR(data,fl-trans_width,fl+trans_width,rp,rs,samplerate,'low');