MATLAB和SciPy爲'buttord'函數提供了不同的結果

我試圖用buttord函數（實際上，我正在移植一個程序，其中這個函數從MATLAB調用到Python）來設計模擬Butterworth過濾器。MATLAB和SciPy爲'buttord'函數提供了不同的結果

我的參數是：

通帶頻率（FP）= 10赫茲，使Wp的= 2 * PI * 10赫茲

阻帶頻率（fs）= 100赫茲，給人Ws的= 2 * PI * 100 Hz

通帶和阻帶損耗/衰減（Rp，Rs）分別爲3和80 dB。

在MATLAB中我使用此代碼：

Wp = 2 * pi * 10 
Ws = 2 * pi * 100 
Rp = 3 
Rs = 80 
[N, Wn] = buttord(Wp, Ws, Rp, Rs, 's')

，讓我N = 5，Wn = 99.581776302。

在SciPy的我試圖做同樣的：

from numpy import pi 
from scipy import signal 
Wp = 2 * pi * 10 
Ws = 2 * pi * 100 
Rp = 3 
Rs = 80 
(N, Wn) = signal.buttord(Wp, Ws, Rp, Rs, analog=True)

，我得到N = 5和Wn = 62.861698649592753。 Wn與MATLAB給出的值不同，並且與Wp非常接近。這裏有什麼問題？我發現this pull request這可能會解釋一些事情：事實證明，MATLAB和SciPy有不同的設計目標（MATLAB試圖針對匹配阻帶頻率進行優化，SciPy試圖針對匹配通帶頻率進行優化）。

我正在使用MATLAB R2013a，Python 3.4.2和SciPy 0.15.0（如果有的話）。

來源

2015-01-20 Renan

Matlab的信號處理工具通常需要頻率相關參數被表示爲相對於該奈奎斯特歸一化頻率。也就是說，Wp = Fp/Fn和Ws = Fs/Fn，其中Fn是奈奎斯特頻率。 – Juderb 2015-01-21 04:58:55

@Juderb我正在使用奈奎斯特頻率無關緊要的模擬濾波器。從SciPy的文檔中，對於模擬濾波器，wp和ws是角頻率（例如rad/s）._ – Renan 2015-01-21 11:18:12

（我也貼了SciPy的郵件列表如下。）

當你設計一個巴特沃斯濾波器buttord，沒有足夠的自由度，以滿足所有的設計約束完全相同。因此是過渡區域的哪一端碰到約束以及哪一端是「過度設計」的選擇。在scipy 0.14.0中做出的改變將該選擇從阻帶邊緣切換到通帶邊緣。

一張照片會說清楚。下面的腳本生成以下圖。（我將Rp從3改爲1.5，-3 dB與Wn的增益一致，這就是爲什麼你的Wn與Wp相同。）使用舊約定或新約定生成的濾波器都滿足設計約束條件。隨着新公約的出現，這種反應恰好與通帶末端的約束相抵觸。

frequency response

import numpy as np 
from scipy.signal import buttord, butter, freqs 
import matplotlib.pyplot as plt 


# Design results for: 
Wp = 2*np.pi*10 
Ws = 2*np.pi*100 
Rp = 1.5  # instead of 3 
Rs = 80 

n_old = 5 
wn_old = 99.581776302787929 

n_new, wn_new = buttord(Wp, Ws, Rp, Rs, analog=True) 

b_old, a_old = butter(n_old, wn_old, analog=True) 
w_old, h_old = freqs(b_old, a_old) 

b_new, a_new = butter(n_new, wn_new, analog=True) 
w_new, h_new = freqs(b_new, a_new) 


db_old = 20*np.log10(np.abs(h_old)) 
db_new = 20*np.log10(np.abs(h_new)) 

plt.semilogx(w_old, db_old, 'b--', label='old') 
plt.axvline(wn_old, color='b', alpha=0.25) 
plt.semilogx(w_new, db_new, 'g', label='new') 
plt.axvline(wn_new, color='g', alpha=0.25) 

plt.axhline(-3, color='k', ls=':', alpha=0.5, label='-3 dB') 

plt.xlim(40, 1000) 
plt.ylim(-100, 5) 

xbounds = plt.xlim() 
ybounds = plt.ylim() 
rect = plt.Rectangle((Wp, ybounds[0]), Ws - Wp, ybounds[1] - ybounds[0], 
        facecolor="#000000", edgecolor='none', alpha=0.1, hatch='//') 
plt.gca().add_patch(rect) 
rect = plt.Rectangle((xbounds[0], -Rp), Wp - xbounds[0], 2*Rp, 
        facecolor="#FF0000", edgecolor='none', alpha=0.25) 
plt.gca().add_patch(rect) 
rect = plt.Rectangle((Ws, ybounds[0]), xbounds[1] - Ws, -Rs - ybounds[0], 
        facecolor="#FF0000", edgecolor='none', alpha=0.25) 
plt.gca().add_patch(rect) 

plt.annotate("Pass", (0.5*(xbounds[0] + Wp), Rp+0.5), ha='center') 
plt.annotate("Stop", (0.5*(Ws + xbounds[1]), -Rs+0.5), ha='center') 
plt.annotate("Don't Care", (0.1*(8*Wp + 2*Ws), -Rs+10), ha='center') 

plt.legend(loc='best') 
plt.xlabel('Frequency [rad/s]') 
plt.ylabel('Gain [dB]') 
plt.show()

來源

2015-01-21 17:15:13

MATLAB和SciPy爲'buttord'函數提供了不同的結果

回答

相關問題