Could anyone help me with my Butterworth filter?

Hey there!

I wanted to use a bandpass filter 1-100Hz (butterworth) and wrote the script below. I'm new to matlab and a little confused, because it appears the filter is still rising at 1Hz en falling from 95-100Hz (it's 0 at 150Hz, the sampling frequency is 300). What I wanted was a flat line at 1 between 1-110 and a rise/fall outside this band quickly to zero. I was wondering if anyone could tell me what I'm doing wrong? The only thing I can think of is using a higher order (which leads me to the question: what order is '' too high '' and what are the negatives of a steep transition band?). I also used FFT to visualize the frequencies to see if it worked (and to practice) but obviously (given the filter) it did not work.

Thanks in advance for any help or explanation!

% script bandpass filter
fs = 300;
fn = 300/2;
fc = [1 110];
[b,a] = butter(4, fc/fn, 'bandpass');

f = 1:1:300;
h=freqz(b,a,f,fs);
figure()
plot(f, abs(h))

data1 = trial {1}(1,:);
SigFil = filtfilt(b,a,data1);
t = time{1}(1,:);
figure()
plot(t,SigFil)

%fft
F = fft(SigFil);
N = length(SigFil);
P2 = abs(F/N);
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = fs*(0:(N/2))/N;
figure()
plot(f,P1)

Answer 1

[Mathieu NOE]

hello

your filter specs and the realized filter are matching

now it seems to me you're looking for a "brick wall" filter characteristics, which means a very high order filter. Unfortunately , often I see the high order filters are unstable when using the butterworth matlab command.

Notice that when you use filtfilt you apply the filter twice , once in forward and one in backward time , so at the end it's like you have applied an order 8 filter , which is already quite a high order.

what do you expect in terms of signal output ?

FYI, your code a bit modified :

clc
clearvars
% script bandpass filter
fs = 300;
fn = 300/2;
fc = [1 110];
[b,a] = butter(4, fc/fn, 'bandpass');
f = logspace(log10(0.25),log10(0.9*fn),300);
h=freqz(b,a,f,fs);
figure()
H_dB = 20*log10(abs(h)); % convert to dB scale

% - 3 dB crossing points 
threshold = -3; % your value here

[fc_low,thr_pos,fc_high,thr_neg]= crossing_V7(H_dB,f,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 

figure(1)
semilogx(f, H_dB,fc_low,thr_pos,'+r',fc_high,thr_neg,'+g','linewidth',2,'markersize',12);grid on
legend('BP filter FRF','fc low','fc high');
xlabel('Freq(Hz)');
ylabel('Gain (dB)');
ylim([-30 3]);

% you can check that fc_low & fc_high are equal to specs (fc)
fc_low
fc_high

% data1 = trial {1}(1,:);
% SigFil = filtfilt(b,a,data1);
% t = time{1}(1,:);
% figure()
% plot(t,SigFil)
% %fft
% F = fft(SigFil);
% N = length(SigFil);
% P2 = abs(F/N);
% P1 = P2(1:N/2+1);
% P1(2:end-1) = 2*P1(2:end-1);
% f = fs*(0:(N/2))/N;
% figure()
% plot(f,P1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%


% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if ~isempty(ind)

    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)

                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end

    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);

    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];

    % extract the negative slope crossing points 
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end

end

[Mathieu NOE]

Hi

problem solved ?

[Mathieu NOE]

hello again

problem solved ?

[Eva]

Hello!

Sorry for the (very) late response, and I forgot to accept your answer! Thank you for your time writing the modified (better) code!!

I was indeed looking for a brick wall filter, what do you mean when you say these are unstable? In terms of signal output, I want to examine PSD after and especially compare the power between 2-100 Hz which are most important to me (if this is what you wanted to know with your question). I thought that a filter would reduce baseline drift and also remove power of frequencies I'm not interested in (and remove the power of a reoccuring large artefact at 130 Hz).

Would you say the order 8 filter is too high? I experimented with using

fn = 300/2;

fc = [1 120];

[b,a] = butter(3, fc/fn, 'bandpass');

and an order 3 filter, but 120 instead of 110, so I am wondering if that matters/makes it more stable?

[Mathieu NOE]

hello Eva

no problem, I know that often people forget to accept an answer ; It's just that I spend quite some time answering here and there , and not everybody thinks about that simple reward.

back to your problem , in fact why the need to do this bandpass filtering on your signal if you can simply do the PSD computation on the raw signal and just keep the PSD data in the frequency range of interest (2-100 Hz) . You can easily put the out of band PSD values to zero or NaN if it's what you want

No need then to do any time domain filering if the purpose is to look at data in the frequency domain

in your comment above , you have just changed a bit the upper corner frequency , but the order is still "low" with 3

now try to generate a butterworth filter with N above 8 and try filter a signal with that. I would not be surprised if the output of the filter is diverging. I have seen that many times.

[Eva]

I understand! It's also a small effort on my part to accept your answer, although that's probably why it's easy to forget as well haha.

You're right, that is much easier and no need for bandpass filtering. Now that I've read your reasoning, it seems so obvious. Thank you! Just out of curiosity - if I had wanted to look at the data in the time-frequency domain, would it be necessary again to do time domain filtering? I feel like the answer is no (because you could also put the out of band tf values to zero) but I'm not sure.

I tried to generate a couple of butterworth filters N above 8 and the signal did change (at some point it even became completely different).

[Mathieu NOE]

hello again

again, if you can do all your work fom the frequency domain data, I don't see the need to do first a time domain filtering. For me it would be needed if you have to stick to the time domain for whatever reason, but it seems not to be the case here

all the best for the future !

[Mathieu NOE]

hello again

just FYI, about designing high order digital filters :

It is best not to use the transfer function implementation of discrete (digital) filters,

because those are frequently unstable. %Use zero-pole-gain and second-order-section implementation instead. Example :

fc = [500 1500];
fs = 1e4;
[z,p,k] = butter(8,fc/(fs/2));
[sos1] = zp2sos(z,p,k);
figure(1)
freqz(sos1,2^16,fs)

●

[Paul]

We also now have zp2ctf, though I'm not really sure what the significant differnence is beween the SOS and CTF forms.