mat实验室脑电带通滤波器

Question

我正在尝试从10分钟长的采样率为500 Hz的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。我已经有过数值精度的问题，它带有如此大的订单的Butterworth过滤器；即使订单低到8，问题是butter为大订单提供了荒谬的b和a值。检查您的b和a向量，您将看到它们包含大约1e21 (!)

解决方案是使用滤波器的零极点表示，而不是系数(b，a)表示。您可以阅读更多关于这个这里的内容。特别地，

通常，您应该使用[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)