matlab與python scipy.signal中的freqs,freqz頻率分析函數，輸出的w，有時候是角頻率，有時候是真實頻率，容易搞混，這里對比一下。
0. 精要總結：
1) freqs:
matlab, 角頻率，rad.s
python, 角頻率 rad/s ,只能是角頻率。
2) freqz
matlab, 形式為 [h,w] = freqz(b,a,n) 角頻率
形式為 [h,f] = freqz(___,n,fs) 時，頻率輸出形式f為Hz形式，fs為采樣頻率
python scipy 中 w,h =freqz(b,a,worN,fs) ， w的單位與輸入fs相同，fs為歸一化角頻率時，w也為角頻率，fs為采樣頻率，單位Hz時，w也為Hz。
3) 角頻率范圍的區別：
freqs中的角頻率是現實中的量，可以很大，比如1000Hz，對應的角頻率為1000*2*pi ; freqz中的角頻率是數字化的，一般使用時是歸一化的，范圍在 0，2*pi之間。
4) 角頻率與Hz頻率轉化：
freqs的w結果要想用Hz,顯示，可以先 w/2/pi 轉化為 Hz 頻率； freqz中的角頻率如果要轉化為具體的頻率， 因為他是歸一化的，用 0~ pi 的范圍代表 0- fs/2 的頻率范圍，可以用 f=（ w/pi）*（fs/2） 轉化為Hz頻率
1. freqs
1.1 matlab中
freqs 是角頻率w的單位 rad/s，想要變成Hz, 顯示時使用 f = w/2/pi
模擬的freqs不存在歸一化。
a = [1 0.4 1];
b = [0.2 0.3 1];
w = logspace(-1,1);
h = freqs(b,a,w);
mag = abs(h);
phase = angle(h);
phasedeg = phase*180/pi;
subplot(2,1,1)
loglog(w,mag)
grid on
xlabel('Frequency (rad/s)')
ylabel('Magnitude')
subplot(2,1,2)
semilogx(w,phasedeg)
grid on
xlabel('Frequency (rad/s)')
ylabel('Phase (degrees)')
1.2 python scipy.signal 中
freqs 輸出的w也是rad/s，也只能是rad/s 角頻率。不過這個不是歸一化的。模擬的freqs不存在歸一化。
w : ndarray
The angular frequencies at which `h` was computed.
b = [1]
a = [0.125 ,1]
# b(0) *s^0
# s = ----------------
# a(0)*s^1 +a(1)*s^0
from scipy.signal import bilinear,freqs,freqz
import matplotlib.pyplot as plt
import numpy as np
# %% python scipy.signal 中 freqs
wf=np.logspace(-1, 4, 1000)
w,h = freqs(b,a,wf)
plt.semilogx(w,20*np.log10(np.abs(h)))
plt.xlabel('rad/s')
如何轉化為Hz顯示，就是x坐標軸 除以 2*pi
b = [1]
a = [0.125 ,1]
# b(0) *s^0
# s = ----------------
# a(0)*s^1 +a(1)*s^0
from scipy.signal import bilinear,freqs,freqz
import matplotlib.pyplot as plt
import numpy as np
# %% python scipy.signal 中 freqs
wf=np.logspace(-1, 4, 1000)
w,h = freqs(b,a,wf)
# plt.semilogx(w,20*np.log10(np.abs(h)))
# plt.xlabel('rad/s')
plt.semilogx(w/2/np.pi,20*np.log10(np.abs(h)))
plt.xlabel('Hz')
2. freqz
2.1 matlab中
1) 函數形式為
[h,w] = freqz(b,a,n)
時，w輸出為角頻率，且歸一化，即最大的角頻率為 pi (對應fs/2,歸一化處理) 。（n為輸出的點的個數）和freqs中
2) 函數形式為
[h,f] = freqz(___,n,fs)
時，頻率輸出形式f為Hz形式，fs為采樣頻率。
b=1;
a=[0.125 1];
fs=2000;
[bz,az] = bilinear(b,a,fs);
[hz0,wz0] = freqz(bz,az,100); % 100是 n，輸出點的個數
fz0= wz0/pi*fs/2; % 將 歸一化的rad/s 轉化為 實際的采樣頻率
plot(fz0,20*log10(abs(hz0)),'o');
xlabel('Hz')
hold on
[hz1,fz1] = freqz(bz,az,100,fs);
plot(fz1,20*log10(abs(hz1)),'r-','linewidth',3);
legend('wz0' , 'fz1')
hold off
結果：
2.2 Python scipy.signal 中
freqz(b,a=1, worN=512, whole=False, plot=None, fs=2*pi, include_nyquist = False,)
Returns
-------
w : ndarray
The frequencies at which `h` was computed, in the same units as `fs`.
By default, `w` is normalized to the range [0, pi) (radians/sample).
w的單位和輸入fs的單位相同，如果fs是用的 rad/s則返回w也是rad/s， 若輸入fs的單位是 Hz,那么輸出的w單位也是Hz。
代碼部分
from scipy.signal import bilinear,freqs,freqz
import matplotlib.pyplot as plt
import numpy as np
# b(0) *s^0
# s = ----------------
# a(0)*s^1 +a(1)*s^0
b = [1]
a = [0.125 ,1]
# %% python scipy.signal 5000中 freqs
wf=np.logspace(-1,4,1000 )
w,h = freqs(b,a,wf)
plt.semilogx(w/2/np.pi,20*np.log10(np.abs(h)),'o',label='freqs')
plt.xlabel('Hz')
plt.ylabel('dB')
fs=5000
bz,az = bilinear(b,a,fs)
worN=np.logspace(-1,4,2000)
idx_end = np.nonzero(worN<=fs/2)[0][-1]
z = freqz(bz,az,worN=worN[0:idx_end],fs=fs)
plt.semilogx(z[0],20*np.log10(z[1]),'-',label='freqz')
plt.legend()
