谈谈MATLAB数字滤波器频域响应
BB: 这几天在做脑电数据的处理工作，遇到很多滤波器上的障碍，在这个函数上也折腾了很久，故记录下一些零碎的感悟。
关于数字滤波器频域响应函数freqz()Abstract: 通常当我们设计好数字滤波器后，会查看其频域响应从而判断该设计是否达到期望指标。频域响应包括其幅度响应和相位响应，可以利用freqz()函数求解。该函数不带返回值时可直接绘制频响特性，接下来我将谈谈该函数带返回值时如何得到一致响应曲线，以及如何利用phasez()函数计算其相位响应。
Note: freqz()、angle()和phasez()简介可看帮助文件MATLAB help center
freqz()函数以官网上给出的一个三阶低通IIR滤波器为例，计算其频率响应。下面是其传递函数：
1. 设计滤波器
根据传递函数得到分母和分析的多项式系数
b0 = 0.05634;
b1 = [1 1];
b2 = [1 -1.0166 1];
a1 = [1 -0.683];
a2 = [1 -1.4461 0.7957];
b = b0*conv(b1,b2); % numerator
a = conv(a1,a2); % denominator
2. 无返回值调用
直接调用freqz(b, a)可以得到频响曲线。更多输入参数设置可以看帮助文件，不同的参数设置可以更改图像的横轴显示方式，如freqz(___,n,fs)，横轴可以显示为对应采样率fs下的实际频率f。
横坐标表示归一化角频率，范围在0 ~ π之间。如果输入参数中带参数’whole’，那么该范围变为0 ~ 2π。需要注意的是：幅度谱纵坐标为分贝，相位谱为角度。
3. 有返回值调用
采用有返回值的调用方式[h,w]=freqz(b, a);其中，h表示滤波器频率响应，为复数形式。如果输入参数中没有设置长度参数n，默认512个点（n的取值范围可以看帮助文档）；w表示角频率，长度和h保持一致，未归一化！我们可以根据h计算频率响应图，为了和无返回值输出结果保持一致，利用abs(h)求其幅度，转化为分贝形式20*log10(abs(h))；利用angle(h)求其相位，此时的角度表示形式为弧度制，还需要转化为角度制，即：angle(h)*180/pi；横坐标以w表示，但需要对其进行归一化，即：w/pi。
hf=abs(h);
hx=angle(h)*180/pi;
subplot(211), plot(w/pi,20*log10(hf)), grid
xlabel('Normalized frequency(*π rad/sample)'), ylabel('Amplitude(dB)')
subplot(212),plot(w/pi,hx), grid
xlabel('Normalized frequency(*π rad/sample)'), ylabel('Phase(degrees)')
可以发现相位响应和上面结果不一样，这是由于angle()函数返回值只能位于±π导致，中间实际上进行了2π（360°）的翻转。
我们也可以利用MATLAB自带的phasez()函数计算滤波器的相位响应曲线，和freqz()函数一样，该函数无返回值时直接输出图像，但是坐标值需要进行转换。调用[phi, w1]=phasez(b, a)，返回值phi以弧度制表示，w为角频率。
[phi, w1]=phasez(b, a)
plot(w/pi, phi*180/pi), grid
xlabel('Normalized frequency(*π rad/sample)'), ylabel('Phase(degrees)')
可以看到，该结果和freqz()函数相位响应保持一致。
%% 3rd-order IIR lowpass filter described by the following polynomial convolutions
b0 = 0.05634;
b1 = [1 1];
b2 = [1 -1.0166 1];
a1 = [1 -0.683];
a2 = [1 -1.4461 0.7957];
b = b0*conv(b1,b2);
a = conv(a1,a2);
%% Data Visualization
% Function (1): freqz
figure(1)
freqz(b, a)
[h,w]=freqz(b, a);
hf=abs(h);
hx=angle(h)*180/pi;
figure(2)
subplot(211), plot(w/pi,20*log10(hf)), grid
xlabel('Normalized frequency(*π rad/sample)'), ylabel('Amplitude(dB)')
subplot(212),plot(w/pi,hx), grid
xlabel('Normalized frequency(*π rad/sample)'), ylabel('Phase(degrees)')
% Function (2): phasez
[phi, w1]=phasez(b, a)
figure(3)
plot(w/pi, phi*180/pi), grid
xlabel('Normalized frequency(*π rad/sample)'), ylabel('Phase(degrees)')
作者：hi强森
