1
$\begingroup$

I want to apply a butterworth bandpass filter to EEG signals and it's not working so I tried with a simple signal:

t=[0:0.01:10];
s3= 9*sin(2*pi*4*t);

Still there is a problem(probably obvious) that i cannot detect. :/ Here is how I design and apply the filter with passband [3 5]Hz:

%Filtering
Fs=100;
Wp1=3;
Wp2=5;
[n,Wn]=buttord([Wp1 Wp2]/(Fs/2), [Wp1-0.5 Wp2+0.5]/(Fs/2), 4, 40); 
[z, p, k] = butter(n, Wn);
[sos,g] = zp2sos(z,p,k);
Xfilt_s3= filtfilt(sos, g, s3);

This is the initial signal and the signal after butterworth application in time and frequency domain. filtered signal in time and frequency

I expected that the signal wouldn't change after the filtering, because its frequency is in the passband. Also, why is this amplitude change occurring? Lastly, since I want to apply it at EEG signals, I would like to know if filtfilt function is appropriate for EEG filtering. Thank you! :)

$\endgroup$
  • 1
    $\begingroup$ Try looking at the transfer function. Some of the filter order routines in Matlab are quirky.Also, butter can directly produce [b a] suitable for filtfilt, you don't need to go from pole zero to sos. Your filter is narrow and filtfilt is equivalent to passing the filter twice. My guess is that your filter is the problem $\endgroup$ – Stanley Pawlukiewicz Jul 30 '17 at 9:05
  • $\begingroup$ Maybe you are looking at the initial transient response? If so, you may have to wait until it settles to view the steady-state response. This time is directly proportional to the filter order. $\endgroup$ – a concerned citizen Jul 31 '17 at 6:06
  • $\begingroup$ I have provided you an answer. It should be solving your problem. $\endgroup$ – Fat32 Jul 31 '17 at 13:02
0
$\begingroup$

You can solve your problem in a number of ways:

First of all you shall better use the a,b coefficients form instead of the z,p,k pole-zero form for designing and applying the filter. Then it can be seen that problem is about the unstability of the designed filter due to very high order with the order returned from the order estimation function.

In your program, your order returns as $n=14$ from the buttord function, with the supplied pass-stop band specifications. However when this order $n$ is used in the butter filter design function with a bandpass type (implicitly specified in the Wn variable which is also returned from buttord function) then the designed filter will have an order of $2n=28$ which is too high and as you can see from its pole-zero plot it's unstable as well:

% S0 - Filter specs
Fs = 100;   % Sampling rate in Hz
Wp1 = 3;    % Pass Band freq1
Wp2 = 5;    % Pass Band freq2

% S1 - Generate the signal (4 Hz sinusoidal)
t=[0:0.01:10];
x= 9*sin(2*pi*4*t);

% S2 - Estimate the order using buttord
[n,Wn]=buttord([Wp1 Wp2]/(Fs/2), [Wp1-0.5 Wp2+0.5]/(Fs/2), 3, 40);

% S3 - Design the filter by the "estimated" order
[b,a] = butter(n,Wn);       % using Wn returned from order estimate

figure,zplane(b,a);         % display the pole-zero plot of the filter
                            % order = 28 , it's untable.
figure,freqz(b,a)           % The frequency plot

See the pole-zero plot (and frequency plot) below: enter image description here enter image description here

The solution is to reduce this order as much as possible. At this point it may not meet the specifications anymore but it will be stable nevertheless. I just divided $n$ by $2$ before calling the butter design function.

[b,a] = butter(n/2,Wn);   
figure,zplane(b,a)          % Now it's stable...(14th order)  
figure,freqz(b,a)           % The frequency plot

See the pole-zero plot (and frequency plot) below: Note that poles are dangerously close to unit circle. Yet still we call it stable. enter image description here enter image description here

Lets apply the stable filter to see if it works as expected:

y1 = filtfilt(b,a,x);

The result is: enter image description here

Note that you don't have to use the filt filt function unless you want zero phase filtering property.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.