Could anyone help me with my Butterworth filter?

조회 수: 9(최근 30일)
Eva 2021년 10월 6일
댓글: Mathieu NOE 2021년 10월 26일
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)

댓글을 달려면 로그인하십시오.

채택된 답변

Mathieu NOE 2021년 10월 7일
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
댓글 수: 6표시숨기기 이전 댓글 수: 5
Mathieu NOE 2021년 10월 26일
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 !

댓글을 달려면 로그인하십시오.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by