【MATLAB数字信号处理】数字滤波器的设计

滤波器指标通常按照其幅度响应指定。例如，低通滤波器$G(z)$的幅度响应$|G(e^{jw})|$通常如图所示，在由$0\le w\le w_p$所定义的通带中，我们要求：

换句话说，幅度以误差$\pm {\sigma }_p$逼近于1。在由$w_s\le |w|\le \pi$定义的阻带中，我们要求：

即幅度以误差${\sigma }_s$接近于零。频率$w_p$和$w_s$分别称为通带边界频率和阻带边界频率，在通带和阻带内的最大容限，${\sigma }_p$和${\sigma }_s$称为波纹。

在大多数应用中，给定的数字滤波器指标如下图所示：

此时，在由$0\le w\le w_p$定义的通带中，幅度的最大值和最小值分别是$1$和$1/\sqrt{1+{ϵ}^2}$，单位为dB的峰通带波纹是

由$w_s\le |w|\le \pi$定义的阻带中，最大波纹由$1/A$表示；单位为dB的最小阻带衰减为

若通带边界频率$F_p$和阻带边界频率$F_s$与数字滤波器的抽样率$F_{\tau }$一起以Hz为单位给定，则以弧度为单位的归一化角边界频率为：

滤波器设计过程的第一步是估计传输函数阶数。对于基于模拟低通滤波器$H_a(s)$的转换的无限冲激响应数字低通滤波器$G(z)$的设计，存在估计滤波器阶数的一个解析公式。对于有限冲激响应低通或高通数字滤波器的设计，可以通过几种设计公式，由下列数字滤波器指标来直接估计最小滤波器长度$N$；归一化通带边界角频率$w_p$；归一化阻带边界角频率$w_s$；峰通带波纹${\sigma }_p$；峰阻带波纹${\sigma }_s$。由Kaiser生成的一个相对简单的逼近式为：

其中$\Delta w=|w_p-w_s|$是过渡带的宽度。上述公式也可用于设计多过渡带有限冲激响应滤波器，在这种情况下，$\Delta w$是所有过渡带的最小宽度。对于具有不等过渡带的多频带滤波器，用上述公式设计的滤波器可能在那些较宽的过渡带中表现出令人无法接受的幅度响应，在这种情况下，这些频带将被处理得比较小，直到得到一个可接受的幅度响应。
​广泛用于设计无限冲激响应滤波器的方法是基于从$s$平面到$z$平面的双线性变换，即

使用上述变换，模拟传输函数$H_a(z)$按照下式被转换成数字传输函数$G(z)$：

通过简单截断前一节给出的理想滤波器的冲激响应系数得到的因果有限冲激响应滤波器，在它们各自的幅度相应中会出现一种摆动行为，称为吉布斯现象。吉布斯现象可用如下方式减弱：用一个两端平滑地降为零的窗口来加窗，或提供一个从通带到阻带平滑的转换。使用渐进窗口会引起旁瓣高度的减小和响应主瓣宽度的增加，得到不连续点处的一个较宽过渡。在所有基于窗口的低通滤波器设计中，截止频率$w_{ϵ}$是通带和阻带边界频率总和的一半。

常用的固定窗函数：
​① 矩形窗

② Barlett窗

③ Hann窗

⑤ Blackman窗

可调节窗函数（这里主要介绍Kaiser窗）：
​Kaiser窗：

式中，$\beta$是可调参数，而$I_o(\mu )$是修正的零阶贝塞尔函数，可用幂级数的形式表示为：

% 这里主要介绍基于窗函数的数字滤波器设计方法
% 载入txt格式的信号幅度数据
data = load('data2.txt');
figure(1);plot(data);grid on
xlabel('time');ylabel('amplitude');title('时域原始信号');
w = 0:0.0114:2*pi;
yf = fft(data);
figure(2);plot(w/pi,abs(yf));grid on
xlabel('frequency');ylabel('amplitude');title('信号的幅度谱');

% 添加随机噪声
scale = 4;
noise = rand(552,1) - 0.5;
noise = scale*noise;
yy = data + noise;
yyf = fft(yy);
figure(3);plot(w/pi,abs(yf),'b');hold on
plot(w/pi,abs(yyf),'r');grid on
xlabel('frequency');ylabel('amplitude');
legend('原始信号','合成信号');
% 将幅度谱变换为分贝
db = 20*log10(abs(yf));
figure(4);plot(db);grid on
xlabel('db');ylabel('amplitude');title('将幅度谱变换为分贝');

%% 矩形窗
wp = 0.2*pi;ws = 0.4*pi;
transWidth = ws - wp;
M1 = 0.92*pi/transWidth;
N1 = ceil(2*M1 + 1);
% N1 = 11
Rwin = rectwin(N1 + 1);
wc1 = (wp + ws)/2;
Wn1 = wc1/pi;
b1 = fir1(N1,Wn1,Rwin);
[h1,omega1] = freqz(b1,1,512);
figure(5);plot(omega1/pi,abs(h1));grid on;
xlabel('\omega/\pi'); ylabel('Magnitude');title('矩形窗的幅度谱')
figure(6);plot(omega1/pi,20*log10(abs(h1)));grid on;
xlabel('\omega/\pi'); ylabel('Gain/dB');title('矩形窗的增益图')
y1 = filtfilt(b1,1,yy);
figure(7);plot(y1,'r');hold on
plot(data);grid on
legend('矩形窗滤波器滤波结果','时域原始信号');

%% Hamming窗
wp = 0.2*pi;ws = 0.4*pi;
transWidth = ws - wp;
M2 = 3.32*pi/transWidth;
N2 = ceil(2*M2 + 1);
% N2 = 35
Hamwin = hamming(N2 + 1);
wc2 = (wp + ws)/2;
Wn2 = wc1/pi;
b2 = fir1(N2,Wn2,Hamwin);
[h2,omega2] = freqz(b2,1,512);
figure(8);plot(omega2/pi,abs(h2));grid on;
xlabel('\omega/\pi'); ylabel('Magnitude');title('Hamming窗的幅度谱')
figure(9);plot(omega2/pi,20*log10(abs(h2)));grid on;
xlabel('\omega/\pi'); ylabel('Gain/dB');title('Hamming的增益图')
y2 = filtfilt(b2,1,yy);
figure(10);plot(y2,'r');hold on
plot(data);grid on
legend('Hamming窗滤波器滤波结果','时域原始信号');

%% Kaiser窗
fpts = [0.2 0.4];
mag = [1,0];
dev = [0.01 0.01];  % 容许误差
[N,Wn,beta,ftype] = kaiserord(fpts,mag,dev);
% N = 23;Wn = 0.3000;beta = 3.3953;ftype = low
Kwin = kaiser(N + 1,beta);
b = fir1(N,Wn,Kwin);
[h,omega] = freqz(b,1,512);
figure(11);plot(omega/pi,abs(h));grid on;
xlabel('\omega/\pi'); ylabel('Magnitude');title('Kaiser窗的幅度谱')
figure(12);plot(omega/pi,20*log10(abs(h)));grid on;
xlabel('\omega/\pi'); ylabel('Gain, dB');title('Kaiser窗的增益图')
y = filtfilt(b,1,yy);
figure(13);plot(y,'r');hold on
plot(data);grid on
legend('Kaiser窗滤波器滤波结果','时域原始信号');

%% 对比三个窗函数的增益图和相位谱
figure(14);
freqz(b1,1);title('矩形窗的增益图和相位谱');
figure(15);
freqz(b2,1);title('Hamming窗的增益图和相位谱');
figure(16);
freqz(b,1);title('Kaiser窗的增益图和相位谱');

%% 原始信号与添加噪声后的合成信号的时域对比
figure(17);plot(data,'b');hold on
plot(yy,'r');grid on
xlabel('time');ylabel('amplitude');
legend('原始信号','合成信号');
figure(18);plot(data,'b');hold on
plot(yy,'r');axis([0 200 -25 25]);grid on
xlabel('time');ylabel('amplitude');
legend('原始信号','合成信号');

figure 1. 原始信号与合成信号（时域局部）

figure 2. 原始信号与合成信号（时域）

figure 3. 原始信号与合成信号（频域）

figure 4. 原始信号的幅度谱

figure 5. 时域原始信号

figure 6. 将幅度谱变换为分贝

figure 7. 矩形窗的幅度谱

figure 8. 矩形窗的增益图

figure 9. 矩形窗的增益图和相位谱

figure 10. 矩形窗滤波结果对比

figure 11. Hamming的幅度谱

figure 12. Hamming窗的增益图

figure 13. Hamming的增益图和相位谱

figure 14. Hamming窗滤波结果对比

figure 15. Kaiser窗的幅度谱

figure 16. Kaiser窗的增益图

figure 17. Kaiser窗的增益图和相位谱

figure 18. Kaiser窗滤波结果对比