matlab - matlab中的EEG带通滤波器

Question

我正在尝试从采样率为 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]);

Answer 1

您的代码(第 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 with zp2sos. 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)