实验室中的EEG带通滤波器 [英] EEG bandpass filter in mat lab

问题描述

我正在尝试从10分钟长的EEG信号(采样率为500Hz)中过滤theta范围(3-8 Hz).这是我的代码.请帮助我了解问题所在.现在，滤波后的信号似乎已损坏.非常感谢！

I'm trying to filter theta range (3-8 Hz) from a 10 min long EEG signal with sampling rate of 500Hz. This is my code. Please help me to understand what's wrong. Right now the filtered signal seems to be ruined. Thank you so much!

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行)给出了等于37的过滤顺序n.我遇到了数值精度的问题>具有如此大订单的巴特沃斯过滤器；即使订单低至8.问题是对于大订单，butter给出了荒谬的b和a值.检查您的b和a向量，您会发现它们包含大约1e21(！)

Your code (line 4) gives a filter order, n, equal to 37. I've had issues of numerical precision with Butterworth filters of such large orders; even with orders as low as 8. The problem is that butter gives absurd b and a values for large orders. Check your b and a vectors, and you'll see they contain values of about 1e21 (!)

解决方案是使用滤波器的零极点表示，而不是系数(b，a)表示.您可以在此处了解更多信息.特别是

The solution is to use the zero-pole representation of the filter, instead of the coefficient (b, a) representation. You can read more about this here. In particular,

通常，您应该使用[z,p,k]语法设计IIR滤波器.要分析或实现您的过滤器，然后可以将[z,p,k]输出与zp2sos一起使用.如果使用[b,a]语法设计过滤器，则可能会遇到数值问题.这些问题是由于舍入错误.对于低至4的过滤器订单，它们可能会发生.

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.

在您的情况下，您可以按照以下步骤进行操作:

In your case, you could proceed along the following lines:

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

这篇关于实验室中的EEG带通滤波器的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！