使用freqz时Octave与SciPy绘图翻转的原因及解决方法

问题背景 #

给定以下滤波器系数：

b=[1.01063287e+00, -1.46490341e+01,  9.94030209e+01, -4.19168764e+02, ...
  1.22949513e+03, -2.66000588e+03,  4.39112431e+03, -5.64225597e+03, ...
  5.70320516e+03, -4.55022454e+03,  2.85602975e+03, -1.39550096e+03, ...
  5.20372994e+02, -1.43160328e+02,  2.74037105e+01, -3.26098385e+00, ...
  1.81735269e-01];
a=[1.00000000e+00, -1.45159238e+01,  9.86464912e+01, -4.16614074e+02, ...
  1.22391361e+03, -2.65216678e+03,  4.38533779e+03, -5.64421414e+03, ...
  5.71487734e+03, -4.56742504e+03,  2.87187255e+03, -1.40575405e+03, ...
  5.25150201e+02, -1.44741759e+02,  2.77584882e+01, -3.30950845e+00, ...
  1.84797453e-01];

在Octave中执行以下代码可得到正常幅频响应图：

pkg load signal
pkg load control

[h(:,1), w] = freqz(flip(b), flip(a),2048);
plot((w/pi),20*log10(abs(h(:,1))));

而Python中执行以下代码时：

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import freqz

h = np.zeros(2048, float)
w = np.zeros(2048, float)

[h[:], w] = freqz(np.flip(b), np.flip(a), 2048)
  
plt.plot((w/np.pi),20*np.log10(h[:]));
plt.show()

两者结果大致一致，但Python生成的图相比Octave呈现类似两次翻转的效果，移除Python代码中的flip命令后结果仍未改变。

原因分析 #

1. 传递函数系数定义差异：
   • Octave的freqz(b,a)默认传递函数为z的负幂次形式：
     $H(z) = \frac{b(1) + b(2)z^{-1} + ... + b(n)z^{-(n-1)}}{a(1) + a(2)z^{-1} + ... + a(m)z^{-(m-1)}}$
     你在Octave中使用flip(b)和flip(a)，相当于将系数反转，传递函数变为$H(z) = \frac{b(n) + b(n-1)z^{-1} + ... + b(1)z^{-(n-1)}}{a(m) + a(m-1)z^{-1} + ... + a(1)z^{{-(m-1)}}$，即对应$H(z}{-1})$的频率响应。
   • Scipy的freqz(b,a)默认传递函数为z的正幂次形式：
     $H(z) = \frac{b[0] + b[1]z + ... + b[n-1]z^{n-1}}{a[0] + a[1]z + ... + a[m-1]z^{m-1}}$
     当你传入翻转后的系数时，得到的是$H(z)$的频率响应，和Octave的$H(z^{-1})$刚好相反，最终导致曲线翻转。
2. 复数处理错误：Python代码中直接对复数h取log10，未取绝对值，虽然结果大致一致，但属于错误操作，可能引发异常。

解决方法 #

方法一：直接使用原始系数（推荐） #

不需要对系数做翻转，直接传入原始b和a，并正确处理复数的幅值：

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import freqz

# 直接传入原始系数，无需翻转
h, w = freqz(b, a, 2048)
# 计算幅值的dB值，和Octave逻辑一致
mag = 20 * np.log10(np.abs(h))
# 绘制幅频响应
plt.plot(w / np.pi, mag)
plt.show()

方法二：使用翻转系数并修正频率响应 #

如果一定要保留翻转系数的操作，可利用实数系数下$H(z^{{-1})$的频率响应等于$H(e}{-j\omega})$共轭的特性，对结果取共轭后再计算幅值：

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import freqz

h, w = freqz(np.flip(b), np.flip(a), 2048)
# 取共轭后计算幅值dB值
mag = 20 * np.log10(np.abs(np.conj(h)))
plt.plot(w / np.pi, mag)
plt.show()

内容的提问来源于stack exchange，提问作者Juha P