如何在没有freqz函数的MATLAB中手动绘制带通Butterworth滤波器的频率响应？

[b, a] = butter(order, [flo fhi]);
filtered_signal = filter(b, a, unfiltered_signal)

#1 楼

我们知道，一般而言，过滤器的传递函数由下式给出：

$$ H（z）= \ dfrac {\ sum_ {k = 0} ^ {M} b_kz ^ {-k}} {\ sum_ {k = 0} ^ {N} a_kz ^ {-k}} $$

现在替换$ z = e ^ {j \ omega} $来评估单元上的传递函数圈子：

$$ H（e ^ {j \ omega}）= \ dfrac {\ sum_ {k = 0} ^ {M} b_ke ^ {-j \ omega k}} {\ sum_ {k = 0} ^ {N} a_ke ^ {-j \ omega k}} $$

因此，这仅是给定$ \ omega $时多项式求值的问题。步骤如下：为频谱的前半部分创建角频率为$ \ omega = [0，\ ldots，\ pi] $的矢量（无需上移至$ 2 \ pi $）并将其保存在w中。
对所有它们都预先计算指数$ e ^ {-j \ omega} $并将其存储在变量ze中。
使用polyval函数进行计算通过调用polyval(b, ze)来计算分子和分母的值，将它们除以并存储在H中。因为我们对振幅感兴趣，所以取结果的绝对值。
通过使用$ H_ {dB} = 20 \ log_ {10} H $转换为dB标度-在这种情况下，$ 1 $是参考值。

将所有内容都放入代码中：

%% Filter definition
a = [1 -0.5 -0.25]; % Some filter with lot's of static gain
b = [1 3 2];

%% My freqz calculation
N = 1024; % Number of points to evaluate at
upp = pi; % Evaluate only up to fs/2
% Create the vector of angular frequencies at one more point.
% After that remove the last element (Nyquist frequency)
w = linspace(0, pi, N+1); 
w(end) = [];
ze = exp(-1j*w); % Pre-compute exponent
H = polyval(b, ze)./polyval(a, ze); % Evaluate transfer function and take the amplitude
Ha = abs(H);
Hdb  = 20*log10(Ha); % Convert to dB scale
wn   = w/pi;
% Plot and set axis limits
xlim = ([0 1]);
plot(wn, Hdb)
grid on

%% MATLAB freqz
figure
freqz(b,a)

freqz的原始输出：



我的脚本输出：



线性比例的快速比较-看起来很棒！

[h_f, w_f] = freqz(b,a);
figure
xlim = ([0 1]);
plot(w, Ha) % mine
grid on
hold on
plot(w_f, abs(h_f), '--r') % MATLAB
legend({'my freqz','MATLAB freqz'})

现在，您可以将其重写为某些函数，并添加一些条件以使其更有用。


另一种方法（以前建议使用的基本特性），即使用滤波器的频率响应是其脉冲响应的傅立叶变换：

$$ H（\ omega）= \ mathcal {F} \ {h（t）\} $$

因此，您必须将$ \ delta（t）$信号馈入系统，计算滤波器的响应并对其进行FFT：

d = [zeros(1,length(w_f)) 1 zeros(1,length(w_f)-1)];
h = filter(b, a, d);
HH = abs(fft(h));
HH = HH(1:length(w_f));

公司mparison这将产生以下结果：

#2 楼

这只是jojek答案的附录，使用双精度数学时，它更通用，也非常好。当精度较低时，当频率响应中的频率非常低（远低于Nyquist）并且滤波器的谐振频率非常低时，就会出现“余弦问题”。

当您计算幅度（平方）频率响应$ | H（e ^ {j \ omega}）| ^ 2 $时，这些复杂的指数将转换为正弦和余弦，但是当进行数学运算时，仅余弦将生存。这应该很清楚，因为幅度是偶函数$ | H（e ^ {-j \ omega}）| = | H（e ^ {j \ omega}）| $ w.r.t.频率，只有余弦是偶函数。

考虑以下触发身份：

$$ \ cos（\ omega）\ = \ 1-2 \ sin ^ 2 \ left （\ frac {\ omega} {2} \ right）$$

现在，看一下等式的右边，所有关于频率的信息都在$ \ sin ^ 2中\ left（\ frac {\ omega} {2} \ right）$项，它从1中减去，并且这个项变得非常小，从$ \ omega \到0 $。如此之小，以至于当将项（或它是负数）加到1时，该项和该项中的频率信息就会丢失。即使是浮点运算，情况也是如此，但双精度浮点运算的问题就更少了。但是我们当中有些将这个频率响应函数投入使用的人可能没有双精度浮点或任何浮点的资源。

所以，我所做的是使用上面的trig标识并消除了所有余弦项，本质上用看起来像$ \ sin ^ 2 \ left（\ frac {\ omega} {2} \ right）$的项以及一些与其他常数组合的常数来代替它们。我将为您展示针对二阶IIR滤波器（又称“二阶”）的答案：

$$ H（z）= \ frac {b_0 + b_1 z ^ {- 1} + b_2 z ^ {-2}} {a_0 + a_1 z ^ {-1} + a_2 z ^ {-2}} $$

它具有复杂的频率响应

$$ H（e ^ {j \ omega}）= \ frac {b_0 + b_1 e ^ {-j \ omega} + b_2 e ^ {-j2 \ omega}} {a_0 + a_1 e ^ {-j \ omega} + a_2 e ^ {-j2 \ omega}} $$

，平方的平方：

$$ \ begin {align}
| H（e ^ {j \ omega}）| ^ 2＆= \ frac {| b_0 + b_1 e ^ {-j \ omega} + b_2 e ^ {-j2 \ omega} | ^ 2} {| a_0 + a_1 e ^ {-j \ omega} + a_2 e ^ {-j2 \ omega} | ^ 2} \\
\\
＆= \ frac {\ big（b_0 + b_1 \ cos（\ omega）+ b_2 \ cos （2 \ omega）\ big）^ 2 + \ big（b_1 \ sin（\ omega）+ b_2 \ sin（2 \ omega）\ big）^ 2} {\ big（a_0 + a_1 \ cos（\ omega）+ a_2 \ cos（2 \ omega）\ big）^ 2 + \ big（a_1 \ sin（\ omega）+ a_2 \ sin（2 \ omega）\ big）^ 2} \\
\\
＆= \ frac {b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1（b_0 + b_2）\ cos（\ omega）+ 2b_0b_2 \ cos（2 \ omega）} {a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1（a_0 + a_2）\ cos（\ omega）+ 2a_0a_2 \ cos（2 \ omega）} \\
\ end {align} $$

所以，一个可以看到幅度频率响应$ | H（e ^ {j \ omega}）| $是偶数对称函数，仅取决于余弦$ \ cos（\ omega）$和$ \ cos（2 \ omega）$。对于非常低的$ \ omega $，这些余弦值非常接近$ 1 $，以单精度定点或浮点数表示，几乎没有剩余的比特可以将这些值与$ 1 $区分开。这就是“余弦问题”。

使用上面的触发身份，您将得到平方的平方：

$$ \ begin {align}
| H（ e ^ {j \ omega}）| ^ 2＆= \ frac {b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1（b_0 + b_2）\ cos（\ omega）+ 2b_0b_2 \ cos（2 \ omega）} {a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1（a_0 + a_2）\ cos（\ omega）+ 2a_0a_2 \ cos（2 \ omega）} \\
\\
＆= \ frac {b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1（b_0 + b_2）\ left（1-2 \ sin ^ 2 \ left（\ tfrac {\ omega} {2} \ right）\ right）+ 2b_0b_2 \ left（1-2 \ sin ^ 2（\ omega）\ right）} {a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1（a_0 + a_2）\ left（1-2 \ sin ^ 2 \ left（ \ tfrac {\ omega} {2} \ right）\ right）+ 2a_0a_2 \ left（1-2 \ sin ^ 2（\ omega）\ right）} \\
\\
＆= \ frac {b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1（b_0 + b_2）\ left（1-2 \ sin ^ 2 \ left（\ tfrac {\ omega} {2} \ right）\ right ）+ 2b_0b_2 \ left（2 \ cos ^ 2（\ omega）-1 \ right）} {a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1（a_0 + a_2）\ left（1-2 \ sin ^ 2 \ left（\ tfrac {\ omega} {2} \ right）\ right）+ 2a_0a_2 \ left（2 \ cos ^ 2（\ omega）-1 \ right）} \\
\\
＆= \ frac {b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1（b_0 + b_2）\ left（1-2 \ sin ^ 2 \ left（\ tfrac {\ omega} {2} \ right）\ right ）+ 2b_0b_2 \ left（2 \ left（1-2 \ sin ^ 2 \ left（\ tfrac {\ omega} {2} \ right）\ right）^ 2-1-right）} {a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1（a_0 + a_2）\ left（1-2 \ sin ^ 2 \ left（\ tfrac {\ omega} {2} \ right）\ right）+ 2a_0a_2 \ left（2 \ left（1 -2 \ sin ^ 2 \ left（\ tfrac {\ omega} {2} \ right）\ right）^ 2-1 \ right）} \\
\\
＆= \ frac {b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1（b_0 + b_2）（1-2-phi）+ 2b_0b_2 \ left（2（1-2 -phi）^ 2-1-\ 1 \ right）} {a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1（a_0 + a_2）（1-2-phi）+ 2a_0a_2 \ left（2（1-2 -phi）^ 2-1-\ right）} \\
\\
＆= \ frac {b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1（b_0 + b_2）（1-2 \ phi）+ 2b_0b_2（1-8 \ phi + 8 \ phi ^ 2）} {a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1（a_0 + a_2）（1 -2 \ phi）+ 2a_0a_2（1-8 \ phi + 8 \ phi ^ 2）} \\
\\
＆= \ frac {b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1b_0 + 2b_1b_2-4（b_1b_0 + b_1b_2）\ phi + 2b_0b_2-16b_0b_2 \ phi + 16b_0b_2 \ phi ^ 2} {a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1a_0 + 2a_1a_2-4（a_1a_0 + a_1a_2）\ phi + 2a_0a_2-16a_0a_2 \ phi + 16a_0a_2 \ phi ^ 2} \\
\\
＆= \ frac {\ big（b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1b_0 + 2b_1b_2 + 2b_0b_2 \ big）-4（b_1b_0 + b_1b_2-4b_0b_2）\ phi + 16b_0b_2 \ phi ^ 2} {\ big（a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1a_0 + 2a_1a_2 + 2a_0a_2 \ big）-4（a_1a_0 + a_1a_2-4a_0a_2）\ phi + 16a_0a_2 \ phi ^ 2} \\
\\
＆= \ frac {\ tfrac14 \ big（b_0 ^ 2 + b_1 ^ 2 + b_2 ^ 2 + 2b_1b_0 + 2b_1b_2 + 2b_0b_2 \ big）-（b_1b_0 + b_1b_2-4b_0b_2）\ phi + 4b_0b_2 \ phi ^ 2} {\ tfrac14 \ big（a_0 ^ 2 + a_1 ^ 2 + a_2 ^ 2 + 2a_1a_0 + 2a_1a_2 + 2a_0a_2 \ big）-（ a_1a_0 + a_1a_2-4a_0a_2）\ phi + 4a_0a_2 \ phi ^ 2} \\
\\
＆= \ frac {\ left（\ frac {b_0 + b_1 + b_2} {2} \ right） ^ 2-\ phi \ big（4b_0b_2（1- \ phi）+ b_1（b_0 + b_2）\ big）} {\ left（\ frac {a_0 + a_1 + a_2} {2} \ right）^ 2-\ phi \ big（4a_0a_2（1- \ phi）+ a_1（a_0 + a_2）\ big）} \\
\ end {align} $$

其中$ \ phi \ triangleq \ sin ^ 2 \ left（\ frac {\ omega} {2} \ right）$

如果您的齿轮打算将其绘制为dB，则显示为

$$ 20 \ log_ {10} | H（e ^ {j \ omega}）| \ = \ 10 \ log_ {10} \ left（\ left（\ tfrac {b_0 + b_1 + b_2} {2} \ right）^ 2-\ phi \ big（4b_0b_2（1- \ phi）+ b_1（b_0 + b_2）\ big）\ right）\\-10 \ log_ {10} \ left（\ left（\ tfrac {a_0 + a_1 + a_2} {2} \ right）^ 2-\ phi \ big（4a_0a_2（1- \ phi）+ a_1（a_0 + a_2）\ big）\ right）$$

，因此您的除法运算将变成减法，但您必须能够计算以某个底数为底的对数。从数字上讲，低频处理的麻烦要比表面上的解决的麻烦少。