matlab - EEG bandpass filter in mat lab

Question

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]);

Answer 1

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 (!)

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,

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)