滤波器设计（2）：经典FIR数字滤波器的设计
引言
以前有一篇IIR滤波器设计的文章。与IIR相对应，这篇文章主要讲FIR滤波器的设计。
以下一段摘自百度百科。
FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，又称为非递归型滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。
详细内容可见FIR滤波器百度词条。
FIR滤波器的优点，一定是稳定的（极点都在0的位置），可以做到线性相位，而IIR的相位是非线性的（摘自郑君里《信号与系统》）。
按奇偶节数、奇偶对称，把实h(n)的FIR分为四类。
从Hg(w)可以看出，有些就不能做低通，有些就不能做高通，等。
其他详细内容，可自行查阅资料。
窗函数法
FIR滤波器设计方法“据说”有很多种。本科上课的时候，就教过一种，窗函数法。
大致意思就是，要设计什么，就用那种对应的理想滤波器去乘以一个窗函数。
假设要设计一个LPF。理想LPF最好理解。但由于xxxx原因，不能物理实现。
上海交大（SJTU）某一年考研题目问的就是“为什么理想LPF物理不可实现”，以及什么是Gibbs（吉布斯）现象。
理想LPF的单位脉冲相应h(n)是无限长的，实际中不可能做到无限长的，所以希望用有限长的h1(n)去逼近h(n)，
实现最小均方误差意义下的最优。
有人证明了，把理想的h(n)截断，就可以得到最小均方误差意义下的最优。
截断的意思就是，加矩形窗w(n)。时域相乘，频域就是拿着H(jw)和窗函数的W(jw)卷积。
后来在有些场合下，发现矩形窗不能满足需求，于是有了三角窗、Hanning窗、Hamming窗、Blackman窗、Kaiser窗等。
一个理解就是拿主瓣宽度和旁瓣衰减互换。鱼和熊掌不可兼得。
e.g. 矩形窗，主瓣窄，旁瓣大，过渡带带宽宽，旁瓣衰减不够大。
总之，需要根据需求，选择窗函数。
用maltab设计FIR线性相位滤波器的步骤大致如下。
1. 根据需求，提出指标。
2. 把指标化成数字的。
3. 按照指标，选窗函数。
4. 按过渡带带宽，定阶数。但是，要注意阶数是否满足频率特性。比如，有些线性相位FIR滤波器就不能实现低通……
5. 调用fir1确定系数。
6. 检验其频率响应。又是一个反复试验，反复修正的过程。
注意，千万小心，matlab里，zplane对FIR滤波器，慎用！慎用！不要轻易对FIR的系数求零极点。
阶数很高的时候，电脑要算死的。
设计实例
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 |
%% 内容：设计FIR滤波器
% 作者：qcy
% 版本：v1.0
% 时间：2016年10月27日22:34:24
clear;
close all;
clc
%%
Fs = 8000; % 采样率
fp = 1000; % Hz 通带截止频率 起伏小于xxx dB
fs = 1500; % Hz 阻带截止频率 衰减大于xxx dB
wp = 2*pi * fp / Fs; % 对应的数字角频率
ws = 2*pi * fs / Fs; % xxx
fc = (fp+fs)/2; % 要设计的FIR的截止频率 -_-! 几个概念容易搞混
% 就LPF来讲，fc按理说可能可以取为 fc = fp + (fs-fp) * theta, theta在0~1之间
w_bw = (ws-wp); % 过渡带带宽
% 如果要求阻带最小的衰减为 50 dB，查表发现hamming窗就可以满足条件
N = ceil(3.3 * 2*pi /w_bw); % 滤波器长度 3.3为Hamming窗过渡段的参数，用来求滤波器长度
win = hamming(N+1);
Fs2 = Fs/2;
Wn = fc/Fs2;
b = fir1(N,Wn,win);
[H,w] = freqz(b,1,1e4);
figure(1);
subplot(211);
stem(b);
title('单位脉冲响应');
grid on;
subplot(212);
plot(w/2/pi*Fs,20*log10(abs(H)+eps));
title('幅频特性');
xlabel('频率（Hz）');
grid on;
figure(2);
subplot(311);
plot(w/2/pi*Fs,20*log10(abs(H)+eps));
title('幅频特性');
grid on;
subplot(312);
plot(w/2/pi*Fs,angle(H));
title('相频特性');
grid on;
subplot(313);
plot(w/2/pi*Fs,unwrap(angle(H)));
grid on;
title('相频特性（解卷绕后）');
xlabel('频率（Hz）');
|
来自CODE的代码片
my_FIR_LPF.m
最后一张图中，相频特性进行了解卷绕（unwrap）。可以看出，相位确实是线性的。