why butterworth bandpass filter changes signal in the passband?

Question

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.

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

Answer 1

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:

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.

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

y1 = filtfilt(b,a,x);

The result is:

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