我正在尝试从采样率为 500Hz 的 10 分钟长 EEG 信号中过滤 theta 范围 (3-8 Hz)。这是我的代码。请帮助我了解问题所在。现在过滤后的信号似乎被破坏了。非常感谢!
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]);
最佳答案
您的代码(第 4 行)给出了一个滤波器阶数 n
,等于 37。我在使用如此大阶数的 Butterworth 滤波器时遇到过数值精度的问题;即使订单低至 8。问题是 butter
为大订单提供荒谬的 b
和 a
值。检查您的 b
和 a
向量,您会看到它们包含大约 1e21 (!)
解决方案是使用滤波器的零极点表示,而不是系数(b
, a
)表示。您可以阅读更多相关信息 here .特别是,
In general, you should use the
[z,p,k]
syntax to design IIR filters. To analyze or implement your filter, you can then use the[z,p,k]
output withzp2sos
. If you design the filter using the[b,a]
syntax, you may encounter numerical problems. These problems are due to round-off errors. They may occur for filter orders as low as 4.
对于您的情况,您可以按照以下几行进行:
[z, p, k] = butter(n,Wn,'bandpass');
[sos,g] = zp2sos(z,p,k);
filt = dfilt.df2sos(sos,g);
fdata = filter(filt,data)
关于matlab - matlab中的EEG带通滤波器,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/23664631/