How to design filter order for high-pass butterworth filter
57 views (last 30 days)
Show older comments
Hi,
I'm new to using filters and was wondering how to choose a filter order. I am currently using a 3rd order and have got the results as below I have seen places that say you can work it out if you have passband edge frequency, stopband edge frequency, passband ripple in dB and stopband attenuation, however I dont understand how you find these for my data. My data set is too large to upload but it is a dynamic movement that looks like as below:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1372669/image.png)
Any help would be appreciated. I have attached my script below if this helps make sense of it.
[Filename,Pathname] = uigetfile('Test 1.txt');
fileID1=fopen([Pathname,Filename]);
formatspec= '%f %f';
C2=textscan(fileID1,formatspec,'HeaderLines',4);
Time2 = C2{1};
Voltage2 = C2{2};
Acceleration2 = (C2{2}/1.995);
Acceleration2 = Acceleration2*9.810;
%*****************************************************************
%define the scanning frequency (Hz) that was used to record the signal
scanfreq= 2048;
n1=length(Acceleration1);
fc= 0.2; %Cut-off Frequency
fs= scanfreq; %sampling Frequency
n_order = 3; %Filter Order
%Ts= delta; %Sampling period
Wn= fc/(fs/2); %Threshold Frequency
%Wn = 0.5;
[b,a] = butter(n_order,Wn,'high');
%****************************************************
%test 1 filtering
Accf1 = detrend(Acceleration1);
Accf1 = filtfilt(b,a,Accf1); %zero-phase digital filtering
Velocity1= cumtrapz(Time1,Accf1); %corrected acceleration Integrated
Displacement1 = cumtrapz(Time1,Velocity1); %Velocity Integrated
[b2,a2] = butter(n_order,Wn,'high');
Velf1=filtfilt(b2,a2,Velocity1); %Filtered Velcoity
DCorrect1= cumtrapz(Time1,Velf1); %Corrected Displacement
DCorrect1= DCorrect1*10^3;
0 Comments
Accepted Answer
Star Strider
on 2 May 2023
Wp= fc/(fs/2); %Threshold Frequency
Ws = 0.95*Wp;
Rp = 1;
Rs = 60;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn,'high')
[sos,g] = zp2sos(z,p,k);
figure
freqz(sos, 2^16, fs)
Velf1=filtfilt(sos,g,Velocity1); %Filtered Velcoity
A transfer function realization can result in an unstable filter. Use a second-order-section representation to avoid that problem.
I would design a bandpass filter for this, because highpass filters tend to select for noise. Having an upper passband lmit can reduce the noise.
.
14 Comments
Star Strider
on 6 May 2023
That definitely helps.
I worked on this for a while. Since the ‘steady state’ regions before and after tha active displacements are not zero (the means of those regions are offset from zero, and integrating those values produces errors in the integration results), so I experimented with setting the beginning and end ‘steady state’ regions before and after the first integrations absolutely to zero. That helps, however it still does not reproduce the LVDT displacement data, although it gets close.
% clc
% clear
JA = readmatrix('Acc.txt');
% JA = JA-0.0154134; %takes it to zero
JA = JA-mean(JA);
JA = detrend(JA);
tA = readmatrix('time.txt');
LVDT = readmatrix('LDVT.txt');
LVDT = LVDT - LVDT(1);
tL = readmatrix('t.txt');
Kistler = readmatrix('AccKis.txt');
% Kistler = Kistler - 0.992151425253108; %takes it to zero
Kistler = Kistler - mean(Kistler);
Kistler = detrend(Kistler);
L = numel(tA);
Fs = 2048;
%JA = detrend(JA);
JA1 = cumtrapz(tA,JA);
JA2 = cumtrapz(tA,JA1);
figure
hold on
plot(tA,JA)
plot(tL,LVDT,'b')
plot(tA,Kistler,'r')
hold off
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTJA = fft(JA.*hann(L), NFFT)/L;
%FTKistler = fft(Kistler.*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
figure
% subplot(3,1,1)
plot(Fv, abs(FTJA(Iv))*2)
grid
title('Fourier transform showing frequency content of JA Accelerometer')
xlabel('frequency (Hz)'); ylabel('Power Spectral Density')
legend('PSD');
hold off
JAsgf = sgolayfilt(JA, 3, 51);
Kistlersgf = sgolayfilt(Kistler, 3, 51);
figure
hold on;
plot(tA, JAsgf)
plot(tA, Kistlersgf)
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG Filtered')
hold off;
FTJAsgf = fft(JAsgf.*hann(L), NFFT)/L;
FTKistlersgf = fft(Kistlersgf.*hann(L), NFFT)/L;
figure
hold on;
plot(Fv, abs(FTJAsgf(Iv))*2)
plot(Fv, abs(FTKistlersgf(Iv))*2)
grid
xlim([0 800])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('JA SG Filtered')
PassbandFreq = 12.5;
xline(PassbandFreq, '-r', 'Filter Passband Frequency')
hold off;
xlim([0 150])
[JAlpf, DF] = lowpass(JAsgf, PassbandFreq, Fs, 'ImpulseResponse','iir'); % 20Hz Passband (Elliptic Filter)
[Kistlerlpf, DFK] = lowpass(Kistlersgf, PassbandFreq, Fs, 'ImpulseResponse','iir');
% DF
figure
hold on;
plot(tA, JAlpf, 'DisplayName','JA')
plot(tA, Kistlerlpf, 'DisplayName','Kistler')
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG & Lowpass Filtered')
hold off
legend('Location','best')
% ylim([-1 1]*0.001)
JAidx(1) = find(JAlpf > 0.3E-3, 1,'first');
JAidx(2) = find(flip(JAlpf) > 0.3E-3, 1,'first');
Kistleridx(1) = find(Kistlerlpf > 0.5E-3,1,'first');
Kistleridx(2) = find(flip(Kistlerlpf) > 0.5E-3,1,'first');
JAlpf(1:JAidx(1)) = 0;
JAlpf(JAidx(2):end) = 0;
Kistlerlpf(1:Kistleridx(1)) = 0;
Kistlerlpf(Kistleridx(2):end) = 0;
JAlpf = JAlpf*10^4; %to scale
Kistlerlpf = Kistlerlpf*10^4;
figure
hold on;
plot(tA(JAidx(1):JAidx(2)), JAlpf(JAidx(1):JAidx(2)), 'DisplayName','JA')
% plot(tA([JAidx(1) JAidx(2)]), JAlpf([JAidx(1) JAidx(2)]), '+r')
plot(tA(Kistleridx(1) : Kistleridx(2)), Kistlerlpf(Kistleridx(1) : Kistleridx(2)), 'DisplayName','Kistler')
% plot(tA([Kistleridx(1) Kistleridx(2)]), Kistlerlpf([Kistleridx(1) Kistleridx(2)]),'xr', 'DisplayName','Kistler')
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG & Lowpass Filtered')
hold off
legend('Location','best')
JV = cumtrapz(tA, JAlpf); % Integrate To Get Velocity
KistlerV = cumtrapz(tA, Kistlerlpf); % Integrate To Get Velocity
JV(1:JAidx(1)) = 0;
JV(JAidx(2):end) = 0;
KistlerV(1:Kistleridx(1)) = 0;
KistlerV(Kistleridx(2):end) = 0;
JD = cumtrapz(tA, JV); % Integrate To Get Displacement
KistlerD = cumtrapz(tA, KistlerV); % Integrate To Get Displacement
figure
hold on;
plot(tA, [JV JD])
plot(tA, [KistlerV KistlerD])
plot(tL,LVDT);
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG & Lowpass Filtered Integrated')
legend('VelocityJA','DisplacementJA','VelocityK','DisplacementK','LVDT','Location','best')
hold off;
ylim([-1 1]*20)
The filters are aboput as effective as they can be, so adjusting their parameters won’t make much difference. I’m not certain what else to suggest at this point.
.
More Answers (0)
See Also
Categories
Find more on Bartlett in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!