Butterworth Bandpass filter implementing
22 次查看(过去 30 天)
显示 更早的评论
Hello everyone,
I tried to filter my EEG Data with the code below. I am not sure if I apply the filter in a wrong way, when I visualize the magnitude spectrum for all channels there is still activity at 50 Hz.
[b,a] = butter(4,[7 35]/(Abtastfrequenz/2),'bandpass');
EEG_filt = filtfilt(b,a,EEG_nurKanaele);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1580056/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1580061/image.png)
0 个评论
回答(1 个)
Star Strider
2024-1-1
There is actually not enough information provided. However in my experience, a generally accepted sampling frequency for EEG is 1 kHz, so using that —
Abtastfrequenz = 1E+3;
[b,a] = butter(4,[7 35]/(Abtastfrequenz/2),'bandpass');
[h,f] = freqz(b, a, 2^16, Abtastfrequenz);
H50 = mag2db(interp1(f, abs(h), 50)) % The Filter Has An Attenuation Of 16.8 dB at 50 Hz
figure
freqz(b, a, 2^16, Abtastfrequenz)
AxMag = subplot(2,1,1);
yl = AxMag.YLim;
hold on
plot(AxMag, [1 1]*50, [yl(1) H50], '-r')
hold off
[z,p,k] = butter(4,[49 51]/(Abtastfrequenz/2),'stop'); % Design A Second Bandstop Filter
[sos,g] = zp2sos(z,p,k); % Second-Order-Section For Stability
[h,f] = freqz(sos, 2^16, Abtastfrequenz);
Hmax = max(abs(h));
H50 = mag2db(interp1(f, abs(h)/Hmax, 50)) % The Filter Has An Attenuation Of 155.9 dB at 50 Hz
figure
freqz(sos, 2^16, Abtastfrequenz)
AxMag = subplot(2,1,1);
yl = AxMag.YLim;
hold on
plot(AxMag, [1 1]*50, [yl(1) H50], '-r')
hold off
My favourite design is an elliptic filter because it is more computationally efficient than other designs, and (in my opinion) performs better. An elliptic design might have enough attentuation at 50 Hz to avoid needing the second bandstop filter, however I am going with your approach (except using the second-order-section implementation in the bandstop filter) in its design and implementation, since using a fourth-order Butterworth design is your choice. You will need to run these filters sequentially (the order is up to you), using the output of the first filter as the input to the second filter. The output of the second filter is the completely filtered signal.
If you have access to the bandpass function, experiment with it with the 'ImpulseResponse','iir' name-value pair to force an elliptic design.
.
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!