matlab實驗中的EEG帶通濾波器

Question

我試圖從採樣率爲500Hz的10分鐘長的EEG信號中過濾θ範圍（3-8Hz）。這是我的代碼。請幫助我瞭解什麼是錯的。目前，過濾的信號似乎已被破壞。非常感謝！

fs=500; 
Wp = [3 8]/(fs/2); Ws = [2.5 8.5]/(fs/2); 
Rp = 3; Rs = 40; 
[n,Wn] = buttord(Wp,Ws,Rp,Rs); 
[b,a] = butter(n,Wn,'bandpass'); 
fdata = filter(b,a,data); 
x=0:ts:((length(data)/fs)-ts); 
f=-fs/2:fs/(length(data)-1):fs/2; 
subplot(2,2,1) 
plot(x,data) 
subplot(2,2,2) 
z1=abs(fftshift(fft(data))); 
plot(f,z1) 
xlim([0 150]); 
subplot(2,2,3) 
plot(x,fdata) 
subplot(2,2,4) 
z=abs(fftshift(fft(fdata))); 
plot(f,z); 
xlim([0 150]); 

Answer 1

您的代碼（第4行）給出了一個濾波器階數，n，等於37。我已經與這樣的大訂單巴特沃斯濾波器的數值精度問題;即使訂單低至8也是如此。問題是butter爲大訂單提供了荒謬的b和a值。檢查b和a載體，你會看到它們包含約1E21值（！）

的解決方案是使用過濾器的零極點表示，代替係數（b，a）表示。你可以閱讀更多關於這個here。特別是，

通常，您應該使用[z,p,k]語法來設計IIR濾波器。要分析或實施您的過濾器，您可以使用輸出與zp2sos。如果使用[b,a]語法設計過濾器，則可能會遇到數值問題。這些問題是由於四捨五入的錯誤。他們可能對濾波器階低至4

在你的情況發生，你可以沿着以下線路進行：

[z, p, k] = butter(n,Wn,'bandpass'); 
[sos,g] = zp2sos(z,p,k); 
filt = dfilt.df2sos(sos,g); 
fdata = filter(filt,data)