matlab中的脑电图带通滤波器

Posted 2023-02-25

【中文标题】matlab中的脑电图带通滤波器

【英文标题】：EEG bandpass filter in mat lab

【发布时间】：2014-05-14 20:36:03

【问题描述】：

我正在尝试从采样率为 500Hz 的 10 分钟长的 EEG 信号中过滤 theta 范围（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]);

【问题讨论】：

张贴图片或示例数据链接，以便我们重现您的问题

【参考方案1】：

您的代码（第 4 行）给出了一个过滤器阶数，n，等于 37。我在使用如此大阶数的 Butterworth 过滤器时遇到了数值精度的问题；即使订单低至 8。问题是 butter 为大订单提供了荒谬的 b 和 a 值。检查您的 b 和 a 向量，您会看到它们包含的值约为 1e21 (!)

解决方案是使用滤波器的零极点表示，而不是系数（b、a）表示。您可以阅读更多关于此here 的信息。特别是，

一般来说，您应该使用[z,p,k] 语法来设计 IIR 滤波器。要分析或实现您的过滤器，您可以使用[z,p,k] 输出和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)

【讨论】：

这是一个很好的解决方案，如果您必须坚持使用原生采样率的 IIR。给定您想要的带通，降低采样率（decimate 或resample）会更好，这样您就可以使用低阶滤波器。或者，改用更稳定的 FIR 滤波器（fir1 或 fir2）。