数字滤波器频率响应分析：从理论到MATLAB实现
数字滤波器频率响应分析：从理论到MATLAB实现
数字滤波器的频率响应分析是信号处理领域的重要内容。本文将介绍如何使用MATLAB中的freqz和freqs函数来计算和可视化数字滤波器的频率响应，包括幅值响应、相位响应等特性。
数字域频率响应
freqz函数使用基于FFT的算法来计算数字滤波器的Z变换频率响应。其基本用法如下：
[h,w] = freqz(b,a,p)
这将返回数字滤波器的包含p个点的复频率响应H(ejω)。其中，向量h包含复频率响应，向量w包含实际频率点。
freqz函数可以接受其他参数，例如采样频率或由任意数量的频率点构成的向量。以下是一个12阶切比雪夫I型滤波器的256点频率响应计算示例：
[b,a] = cheby1(12,0.5,200/500);
[h,f] = freqz(b,a,256,1000);
在这个例子中，采样频率fs被指定为1000Hz。
归一化频率
该工具箱使用单位频率是奈奎斯特频率的约定，定义为采样频率的一半。所有基本滤波器设计函数的截止频率参数均用奈奎斯特频率进行归一化处理。例如，对于采样频率为1000Hz的系统，300Hz等于300/500=0.6。
要将归一化频率转换为围绕单位圆的角频率，请乘以π。要将归一化频转换回赫兹，请乘以采样频率的一半。
幅值和相位响应绘制
如果不带输出参数调用freqz，它会同时绘制幅值对频率的图和相位对频率的图。例如，截止频率为400Hz、基于2000Hz的采样频率的九阶巴特沃斯低通滤波器的频率响应绘制如下：
[b,a] = butter(9,400/1000);
freqz(b,a,256,2000)
freqz也可以接受由任意数量的频率点构成的向量，以用于频率响应计算。例如：
w = linspace(0,pi);
h = freqz(b,a,w);
这将计算由向量b和a定义的滤波器在w的频率点处的复频率响应。频率点可以是0到2π范围内的值。要指定从零到采样频率的频率向量，请在参数列表中同时包括频率向量和采样频率值。
示例：传递函数的频率响应
计算并显示由以下传递函数描述的三阶IIR低通滤波器的幅值响应：
将分子和分母表示为多项式卷积。求分布在整个单位圆上的2001个点上的频率响应。
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);
[h,w] = freqz(b,a,'whole',2001);
绘制以分贝表示的幅值响应：
plot(w/pi,20*log10(abs(h)))
ax = gca;
ax.YLim = [-100 20];
ax.XTick = 0:.5:2;
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')
示例：FIR带通滤波器的频率响应
设计一个FIR带通滤波器，通带在0.35π和0.8π弧度/采样点之间，波纹为3dB。第一个阻带是从0到0.1π弧度/采样点，衰减为40dB。第二个阻带是从0.9π弧度/采样点到奈奎斯特频率，衰减为30dB。计算频率响应。同时以线性单位和分贝绘制其幅值。突出显示通带。
sf1 = 0.1;
pf1 = 0.35;
pf2 = 0.8;
sf2 = 0.9;
pb = linspace(pf1,pf2,1e3)*pi;
bp = designfilt('bandpassfir', ...
'StopbandAttenuation1',40, 'StopbandFrequency1',sf1,...
'PassbandFrequency1',pf1,'PassbandRipple',3,'PassbandFrequency2',pf2, ...
'StopbandFrequency2',sf2,'StopbandAttenuation2',30);
[h,w] = freqz(bp,1024);
hpb = freqz(bp,pb);
subplot(2,1,1)
plot(w/pi,abs(h),pb/pi,abs(hpb),'.-')
axis([0 1 -1 2])
legend('Response','Passband','Location','South')
ylabel('Magnitude')
subplot(2,1,2)
plot(w/pi,db(h),pb/pi,db(hpb),'.-')
axis([0 1 -60 10])
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Magnitude (dB)')
示例：高通滤波器的幅值响应
设计一个三阶高通巴特沃斯滤波器，它具有0.5π弧度/采样点的归一化3-dB频率。计算它的频率响应。用分贝表示幅值响应，并对其绘图。
[b,a] = butter(3,0.5,'high');
[h,w] = freqz(b,a);
dB = mag2db(abs(h));
plot(w/pi,dB)
xlabel('\omega / \pi')
ylabel('Magnitude (dB)')
ylim([-82 5])
使用fvtool重复计算：
fvtool(b,a)
模拟域频率响应
freqs函数计算由两个输入系数向量b和a定义的模拟滤波器的频率响应。其运算类似于freqz的运算；您可以指定要使用的频率点数量，提供由任意数量的频率点构成的向量，并绘制滤波器的幅值和相位响应。
示例：模拟IIR低通滤波器的比较
尝试此示例Copy Code Copy Command
设计截止频率为2GHz的五阶模拟巴特沃斯低通滤波器。乘以2π以将频率转换为弧度/秒。计算滤波器在4096个点上的频率响应。
n = 5;
fc = 2e9;
[zb,pb,kb] = butter(n,2*pi*fc,"s");
[bb,ab] = zp2tf(zb,pb,kb);
[hb,wb] = freqs(bb,ab,4096);
设计一个具有相同边缘频率和3dB通带波纹的五阶切比雪夫I型滤波器。计算它的频率响应。
[z1,p1,k1] = cheby1(n,3,2*pi*fc,"s");
[b1,a1] = zp2tf(z1,p1,k1);
[h1,w1] = freqs(b1,a1,4096);
设计一个具有相同边缘频率和30dB阻带衰减的5阶切比雪夫II型滤波器。计算它的频率响应。
[z2,p2,k2] = cheby2(n,30,2*pi*fc,"s");
[b2,a2] = zp2tf(z2,p2,k2);
[h2,w2] = freqs(b2,a2,4096);
设计一个具有相同边缘频率和3dB通带波纹、30dB阻带衰减的五阶椭圆滤波器。计算它的频率响应。
[ze,pe,ke] = ellip(n,3,30,2*pi*fc,"s");
[be,ae] = zp2tf(ze,pe,ke);
[he,we] = freqs(be,ae,4096);
设计一个具有相同边缘频率的5阶贝塞尔滤波器。计算它的频率响应。
[zf,pf,kf] = besself(n,2*pi*fc);
[bf,af] = zp2tf(zf,pf,kf);
[hf,wf] = freqs(bf,af,4096);
对衰减绘图，以分贝为单位。以千兆赫为单位表示频率。比较滤波器。
plot([wb w1 w2 we wf]/(2e9*pi), ...
mag2db(abs([hb h1 h2 he hf])))
axis([0 5 -45 5])
grid
xlabel("Frequency (GHz)")
ylabel("Attenuation (dB)")
legend(["butter" "cheby1" "cheby2" "ellip" "besself"])
巴特沃斯和切比雪夫II型滤波器具有平坦的通带和宽过渡带。切比雪夫I型和椭圆滤波器转降更快，但有通带波纹。切比雪夫II型设计函数的频率输入设置阻带的起点，而不是通带的终点。贝塞尔滤波器沿通带具有大致恒定的群延迟。