How to design filter order for high-pass butterworth filter

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:

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;

Answer 1

[Star Strider]

To determine the optimal order, use the buttord function.

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.

.

[Emily Keys]

Hi, Thanks so much for this. Also Would it be possible for you to explain why you use the values you do for Ws, Rp and Rs if you wouldnt mind? If you have any websites/ reports that explain filter design I would appreciate it so much also.

Also one last thing, are you suggesting just putting a filter between the velocity and not for the acceleration too? I havent got it to work just yet but thanks again for your response and any other suggestions would be greatly appreciated

scanfreq= 2048;

n=length(Acceleration);

%n=length(y);

fc= 0.2; %Cut-off Frequency

fs= scanfreq; %sampling Frequency

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)

[b,a] = butter(n,Wp,'high');

%***************************************************

%test 2 filtering

Accf = Acceleration;

Velocity = cumtrapz(Time,Accf); %corrected acceleration Integrated

Displacement = cumtrapz(Time,Velocity); %Velocity Integrate

Velocity=filtfilt(sos,g,Velocity);

DCorrect= cumtrapz(Time,Velocity); %Corrected Displacement

DCorrect= (DCorrect*10^3);

[Star Strider]

The values I use for those parameters are simply from my experience. You can use any values you want. The passband is apparently 0.2 Hz (or whatever frequency unit you’re using) and infinity. I chose a stopband (‘Ws’) frequency that was close to it because it is itself extremely low. With frequencies in this range, an elliptic filter woudl likely be a better option. The design approach would be essentially the same (the elliptic filter design functions have slightly different arguments).

I’m not certain what your data are, hwoever you need to be sure that the sampling intervals are all the same. The easiest way to do that is to calculate:

TsMean - mean(diff(t))

TsStd = std(diff(t))

I am assuming here that ‘t’ is your sampling time vector.

If the standard deviation is on the order of or less of the mean, the sampling frequency is likely constant enough. If it is much higher than that, you likely need to use the resample function to provide a regular sampling frequency, because essentially all signal processing (the exception being nufft) requires it. Then do the filtering.

It would likely be appropriate to calculate the Fourier transform (fft or pspectrum) of your signal before processing it. That way, you understand what its frequency components are and can design the filters appropriately.

Also, you might get better results if you subtract the mean of the signal from the enitire signal before doing the Fourier transform or any filtering, then filter it. You can add the mean of the signal back later if necessary.

.

[Emily Keys]

Hi, Thanks for messaging back. I have now started to try apply an elliptic filter to another set of data which involved a quasi-static response and am having some issues fitting it correctly. I wish I could send the data set but unfortunately it is too large but here is the fft and the results I have got so far. As you can see theres a lot more high frequency noise so I have used a bandpass filter, is it just a matter of playing around with the filter parameters until it works? Any advice would be appreciated greatly yet again and I hope I am not asking too much.

Kistler = detrend(Kistler(4368:32582,9)); %shortening data and detrending
JA = detrend(JA(4368:32582,9));
tA = tA(4368:32582,9);

for i = 1 : 1

JAmean(:,i) = mean(JA(:,i));
JA(:,i) = JA(:,i) - JAmean(:,i);  %removing mean from signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% filter stuff
scanfreq = fsA;

n2=length(JA(:,i));


fc= 0.1; %Cut-off Frequency
f2 = 1; 
fs= scanfreq; %sampling Frequency
Wp = [fc f2]/(fs/2); %Threshold Frequency
Ws = [0.09 1.01]/(fs/2); %
Rp = 1;
Rs = 60;


[n,Wn] = ellipord(Wp,Ws,Rp,Rs);

[z,p,k] = ellip(n,Rp,Rs,Wn,'bandpass');

[sos,g] = zp2sos(z,p,k);

freqz(sos, 2, fs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%filter and integrated for JA

AccfJA(:,i) = filtfilt(sos,g,JA(:,i)); %zero-phase digital filtering
VelocityJA(:,i)= cumtrapz(tA(:,i),AccfJA(:,i)); %Velocity Integrated
%Displacement = cumtrapz(tA(:,i),Velocity(:,i));

VelfJA(:,i)=filtfilt(sos,g,VelocityJA(:,i)); %Filtered Velcoity

DCorJA(:,i)= cumtrapz(tA(:,i),VelfJA(:,i)); %Corrected Displacement
DCorJA(:,i)= DCorJA(:,i)*10^4;

%[Pyy1,f1] = FourierTran(scanfreq,JA(:,i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%non-filter and integrated for JA
VelJA(:,i)= cumtrapz(tA(:,i),JA(:,i));
DispJA(:,i) = cumtrapz(tA(:,i),VelJA(:,i));
DispJA(:,i) = DispJA(:,i)*10^4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%filter and integrated for Kistler

AccfK(:,i) = filtfilt(sos,g,Kistler(:,i)); %zero-phase digital filtering
VelocityK(:,i)= cumtrapz(tA(:,i),AccfK(:,i)); %Velocity Integrated
%Displacement = cumtrapz(tA(:,i),Velocity(:,i));

VelfK(:,i)=filtfilt(sos,g,VelocityK(:,i)); %Filtered Velcoity

DCorK(:,i)= cumtrapz(tA(:,i),VelfK(:,i)); %Corrected Displacement
DCorK(:,i)= DCorK(:,i)*10^4;

%[Pyy2,f2] = FourierTran(scanfreq,Kistler(:,i));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%non-filter and integrated for Kistler
VelK(:,i)= cumtrapz(tA(:,i),Kistler(:,i));
DispK(:,i) = cumtrapz(tA(:,i),VelK(:,i));
DispK(:,i) = DispK(:,i)*10^4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plotting graphs


figure();clf; 
    subplot(3,1,1)
    plot(tA(:,i),DCorJA(:,i));hold on
    plot(tA(:,i),DCorK(:,i));
    plot((tL-2),LDVT(:,i));hold on
    %xlim([0 20]);
   % legend('JA','Kistler','LVDT')
    xlabel('Time (s)');   ylabel('Displacement (mm)');
    hold off

[Star Strider]

It’s difficult for me to understand those data. There appear to be about 2 cycles over 15 seconds and it seems to be sampled at about 2400 Hz. There seems to be a peak at about 20 Hz and what may be mains frequency contamination at about 50-60 Hz, and significant high-frequency noise (that I doubt is actual data unless there is some sort of 500 Hz mechanical interference). I have no idea what the actual data are that you want the filter to pass (and it seems that I will never have the actual data to work with, although if you use the resample function to downsample it and reduce its size, it might be possible to upload it).

The PSD plot is helpful, however without being able to explore the data myself, as well as knowing what you want to get from it, I cannot suggest a specific approach.

The data have broadband noise, and the easiest way to deal with that is to use the Savitzky-Golay filter (sgolayfilt function) to filter it out. I usually use a 3-order polynomial, and like all signal processing tasks, it will be necessary to experiment with the ‘framelen’ argument to get the best result.

Design the filter to get the data from the signal that you want. Experiment until you get the result you want. Unfortunately I can’t be more specific than that, without knowing more about what the data are and what you want from them.

[Emily Keys]

Hi, Im sorry about not being able to upload it, but I think I have found a way to cut the data and attach it now which is great. I should have been clearer about what im trying to achieve. I am basically double integrating accelerometer data to get the displacement of a weight going over a model bridge to model a quasi-static response. I am trying to get the data to look like the LVDT/ accurate data that was also measured from the experiment but have not found a method yet that works to an acceptable level of error.

I have used a Savitzy-Golday filter on this data which will show in the code below although sorry as the code is very messy. It worked a well hoewever the post-load became quite odd as can be seen below, thanks again for youre response also:

for test = 1 : 9 %matlab code was for 9 sets of data but im focusing on 1
    
    [~,a] = min(abs(tAcc-AccStart(test)));
    [~,b] = min(abs(tAcc-AccEnd(test)));

    JA(:,test) = JAWhole(a:b);
    Kistler(:,test) = KistlerWhole(a:b);
    
    JA(:,test) = detrend(JA(:,test));
    Kistler(:,test) = detrend(Kistler(:,test));

    tA(:,test) = tAcc(a:b);

    tA(:,test) =  tA(:,test) -( tA(1,test) );
    
end

%Kistler = Kistler(4368:32582,9);
%JA = JA(4368:32582,9);
%tA = tA(4368:32582,9);
Kistler = Kistler(:,9);
JA = JA(:,9);
tA = tA(:,9);


clear KistlerWhole JAWhole tAcc AccStart AccEnd test a b
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;


hFig1 = figure;
%tA = tA(:,9);
%JA = JA(:,9);
plot(tA, JA, 'b.-', 'MarkerSize', 1);
grid on;
hold on;
fontSize = 20;
xlabel('Time', 'FontSize', fontSize);
ylabel('Acceleration', 'FontSize', fontSize);
title('Original Signal', 'FontSize', fontSize);
hFig1.WindowState = 'maximized'; % Maximize the figure window.
% Draw a line at y=0
yline(0, 'LineWidth', 2);

% A moving trend is influenced by the huge outliers, so get rid of those first.
% Find outliers
outlierIndexes = isoutlier(JA);
plot(tA(outlierIndexes), JA(outlierIndexes), 'ro', 'MarkerSize', 15);
% Extract the good data.
tGood = tA(~outlierIndexes);
accelGood = JA(~outlierIndexes);
% plot(t(~outlierIndexes), accel(~outlierIndexes), 'mo', 'MarkerSize', 10); % Plot circles around the good data.
% Do a Savitzky-Golay filter (moving quadratic).
windowWidth = 51; % Smaller for tighter following of original data, bigger for smoother curve.
smoothedy = sgolayfilt(accelGood, 2, windowWidth);
hold on;
plot(tGood, smoothedy, 'r-', 'LineWidth', 2);
%legend('Original Signal', 'X axis', 'Outliers', 'Smoothed Signal');

% fill in the missing points.
smoothedy = interp1(tGood, smoothedy, tA);
% Now subtract the smoothed signal to get the variation
signal = JA - smoothedy;
% Plot it.
hFig2 = figure;
plot(tA, signal, 'b.-', 'MarkerSize', 9);
grid on;
hold on;
title('Corrected Signal', 'FontSize', fontSize);
xlabel('Time', 'FontSize', fontSize);
ylabel('Acceleration', 'FontSize', fontSize);
hFig2.Units = 'normalized';
hFig2.Position = [.2, .2, .5, .5]; % Size the figure window.
% Draw a line at y=0
yline(0, 'LineWidth', 2);


VelCor = cumtrapz(tA,signal);
DispCorr = cumtrapz(tA,VelCor);
DispCorr = DispCorr*10^4;

figure(); clf; hold on
subplot(1,1,1);
plot(tA, DispCorr, 'b');hold on
plot(tL,LDVT(:,9));hold off
grid on;
hold on;
title('Corrected Displacement', 'FontSize', fontSize);
xlabel('Time', 'FontSize', fontSize);
ylabel('Displacement', 'FontSize', fontSize);
hFig2.Units = 'normalized';
hFig2.Position = [.2, .2, .5, .5]; % Size the figure window.
% Draw a line at y=0

[Star Strider]

My pleasure!

I took a slightly different approach, since I do not completely understand your approach. I am not certain how the result is supposed to look (structural engineering is far from my areas of expertise), however it should be relatively straightforward to follow what I did here.

I would not remove the outliers, first because that distorts the signal and second because there is only one significant outlier and that can be dealt with by the Savitzky-Golay filter. (In any event, using filloutliers to interpolate them is preferable to removing them.)

My approach to the initial analysis —

JA = readmatrix('JA.txt');
tA = readmatrix('tA.txt');
L = numel(tA);
Fs = 2048;

figure
plot(tA, JA)
grid
xlabel('Time')
ylabel('Amplitude')
title('JA')

figure
histfit(JA, 500)
grid
xlim([-1 1]*0.15)
title('Check Distribution')

Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTJA = fft(JA.*hann(L), NFFT)/L;
Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTJA(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('JA Unfiltered')

JAsgf = sgolayfilt(JA, 3, 51);

figure
plot(tA, JAsgf)
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG Filtered')

FTJAsgf = fft(JAsgf.*hann(L), NFFT)/L;

figure
plot(Fv, abs(FTJAsgf(Iv))*2)
grid
xlim([0 800])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('JA SG Filtered')
xline(20, '-r', 'Filter Passband Frequency')

[JAlpf, DF] = lowpass(JAsgf, 20, Fs, 'ImpulseResponse','iir');      % 20Hz Passband (Elliptic Filter)
% DF

figure
plot(tA, JAlpf)
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG & Lowpass Filtered')

JV = cumtrapz(tA, JAlpf);                                           % Integrate To Get Velocity
JD = cumtrapz(tA, JV);                                              % Integrate To Get Displacement

figure
plot(tA, [JV JD])
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG & Lowpass Filtered Integrated')
legend('Velocity','Displacement', 'Location','best')

If my analysis is not appropriate or correct, please point out any errors in it to educate me in the correct analytic procedures for these data.

.

[Emily Keys]

Hi, Thanks so much for showing these results, incredible!! I would say it looks really good and I will attached the LVDT data so you can see what it is meant to look like, the peak is meant to be around -24 mm. I would say it looks like it needs to be scaled by a power of 10 which is not a big deal at all but also if you would know any ways to correct the pre load (before about 7ish seconds) and post load (at about 12 ish seconds onwards) as they should just be flat in theoery. But its still amazing work and thank you so much again.

[Star Strider]

Having the LVDT data helps somewhat. The problem is that the LVDT data are sampled at 8.2 Hz, while the accelermoeter signals are sampled at 2048 Hz. This means that there is a significant amount of missing information in the LVDT signal. It would be interesting to know what those data would be if the two signals were sampled at 2048 Hz.

The best I can do is to design a compensator filter using the transfer function defined as where represents the Fourier transform, after resampling the LVDT signal to equal the JA signal (this is not generally appropriate because it creates data where no data previously existed, however it is unavoidable in this instance), and equalising their lengths (this means truncating the JA signal, however since it is essentially zero in that region, no information is lost). I used the original ‘JA’ signal for this analysis, not the filtered versions.

To estimate the filter, I used the Signal Processing Toolbox invfreqz function instead of the System Identification (SI) Toolbox functions because I know you have the Signal Processing Toolbox. The SI functions might give a better result, however considering the lost information in the LVDT signal, that probably doesn’t matter that much. The filter is stable, however it is not minimum-phase.

The compensator filter with the ‘JA’ signal produces the correct waveform, however not the correct amplitude or time-domain characteristics (it produces a significantly longer result than it should, and varying the filter order does not change that). Downsampling (using resample) the 2048 Hz accelerometer signal to the 8.2 Hz LVDT signal might solve this, however that will be at the expense of discarding a significant amount of information, even though no new information would be created. I have not doone that analysis with these data.

This is not the high-quality result that I want, however given the constraints, it’s likely the best that can be done with these data. I encourage you to experiment with it, specifically to design a filter for the ‘JA’ signal the can produce something similar to the ‘LVDT’ signal.

This analysis should get you started —

JA = readmatrix('JA.txt');
tA = readmatrix('tA.txt');
LA = numel(tA);
FsA = 2048;


% figure
% plot(tA, JA)
% grid
% xlabel('Time')
% ylabel('Amplitude')
% title('JA')

% figure
% histfit(JA, 500)
% grid
% xlim([-1 1]*0.15)
% title('Check Distribution')

FnA = FsA/2;
NFFTA = 2^nextpow2(LA);
FTJA = fft(JA.*hann(LA), NFFTA)/LA;
FvA = linspace(0, 1, NFFTA/2+1)*FnA;
IvA = 1:numel(FvA);

figure
plot(FvA, abs(FTJA(IvA))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('JA Unfiltered')

JAsgf = sgolayfilt(JA, 3, 51);

figure
plot(tA, JAsgf)
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG Filtered')

FTJAsgf = fft(JAsgf.*hann(LA), NFFTA)/LA;

figure
plot(FvA, abs(FTJAsgf(IvA))*2)
grid
xlim([0 800])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('JA SG Filtered')
xline(0.8, '-r', 'Filter Passband Frequency')

[JAlpf, DF] = lowpass(JAsgf, 0.8, FsA, 'ImpulseResponse','iir');      % 20Hz Passband (Elliptic Filter)
% DF

figure
plot(tA, JAlpf)
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG & Lowpass Filtered')

JV = cumtrapz(tA, JAlpf);                                           % Integrate To Get Velocity
JD = cumtrapz(tA, JV);                                              % Integrate To Get Displacement

figure
plot(tA, [JV JD])
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG & Lowpass Filtered Integrated')
legend('Velocity','Displacement', 'Location','best')

LVDT = readmatrix('LVDT.txt');
tL = readmatrix('tL.txt').';
tL = tL - tL(1);                                                        % Start At Zero
% format long
% ts = [mean(diff(tL)); std(diff(tL)); std(diff(tL))/mean(diff(tL))]

figure
plot(tL, LVDT)
grid
xlabel('tL')
ylabel('LVDT')

LL = numel(tL);
FsL = 1/mean(diff(tL))

FsL = 8.2000

FnL = FsL/2;
NFFTL = 2^(nextpow2(LL)+6);
FTLVDT = fft((LVDT-mean(LVDT)).*hann(LL), NFFTL)/LL;
FvL = linspace(0, 1, NFFTL/2+1)*FnL;
IvL = 1:numel(FvL);

figure
plot(FvL, abs(FTLVDT(IvL))*2)
grid
xlim([0 1])
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('LVDT Unfiltered')
xline(0.8, '-r', '800 mHz')

[LVDTr,tLr] = resample(LVDT,tL,FsA,'linear');               % Upsample 'LVDT' To 2048 Hz.

% Sz = [size(LVDTr); size(JA)]

JAv = JA(1:numel(LVDTr));
tAv = tLr(1:numel(tLr));

% JAv = JAv * rms(LVDT)/rms(JAv);

figure
plot(tAv, JAv)
grid
xlabel('Time')
ylabel('Amplitude')
title('JAv')

LAv = numel(JAv);
FnAv = FsA/2;
NFFTAv = 2^nextpow2(LAv);
FTJAv = fft(JAv.*hann(LAv), NFFTAv)/LAv;
FvAv = linspace(0, 1, NFFTAv/2+1)*FnAv;
IvAv = 1:numel(FvAv);

figure
plot(FvAv, abs(FTJAv(IvAv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('JAv Unfiltered')

figure
plot(tLr, LVDTr)
grid
xlabel('tL')
ylabel('Amplitude')
title('LVDT Resampled')

LLr = numel(tLr);
FsLr = 1/mean(diff(tLr));
FnLr = FsLr/2;
NFFTLr = 2^nextpow2(LLr);
FTLVDTr = fft((LVDTr-mean(LVDTr)).*hann(LLr), NFFTLr)/LLr;
FvLr = linspace(0, 1, NFFTLr/2+1)*FnLr;
IvLr = 1:numel(FvLr);

Sz = [size(FTLVDTr); size(FTJAv)]

Sz = 2×2

65536 1 65536 1

TFq = FTLVDTr ./ FTJAv;

figure
plot(FvAv, abs(TFq(IvAv)))
grid
xlim([0 1])
xlabel('Frequency (Hz)')
ylabel('Magnitude)')
title('Transfer Function: LVDT/JA')

figure
plot(FvAv, imag(TFq(IvAv)))
grid
xlim([0 5])
ylim([-1 1]*1500)
xlabel('Frequency (Hz)')
ylabel('Magnitude)')
title('Transfer Function: \Im{(LVDT/JA)}')

[b,a] = invfreqz(TFq(IvAv), FvAv*pi/FvAv(end), 5, 5);                       % Design Compensator Filter

figure
freqz(b,a,2^16,FsA)

figure
zplane(b,a)
grid

JAv_filt = filtfilt(b,a,JAv);

figure
plot(tAv, JAv_filt)
grid
xlabel('Time')
ylabel('Amplitude')
title('Filtered ‘JAv’ Signal')

If it is possible to get the LVDT signal to be sampled at 2048 Hz, it would be helpful to have that signal. There might be fewer problems designing the filter than there are with the currently available signals, and alternatively, the two signals might be close enough that only denoising is necessary.

.

[Emily Keys]

Hi, thank you so much again for working on this. Sorry I should have been clear that I don't want to use the LVDT data to aid in the accelerometer data since I am trying to get the data as accurately as possible by just using the accelerometer and filters and baseline correction method. Saying that you have given a great lot to work with so thank you so much for helping!

[Star Strider]

As always, my pleasure!

If you only want to use the accerometer data, experiment with the lowpass cutoff frequencies to see what gives the best results. I would use the signal processing code in my earlier Comment (not all of it is necessary, since some of it like the histfit call was simply my exploring your data) and experiment with this:

[JAlpf, DF] = lowpass(JAsgf, 20, Fs, 'ImpulseResponse','iir'); % 20Hz Passband (Elliptic Filter)

or design your own elliptic filter to get the most appropriate result. Your accelerometer data are high-quality (I checked for that as well), so it should not be necessary to do any processing on it other than what we’ve done here.

I will of course continue to help as necessary if you have any further questions.

[Emily Keys]

Hi again, I was just wondering if you would be able to look over this other data I am working on. i am trying to again apply your method of the savizty-golay filter and lowpass cutoff frequency but cant seem to get it right when tweaking it. I have attached the code and a photo of what it keeps looking like. I have also added a line to detrend the acceleration but it doesnt seem to change it. Any advice/help would be massively appreciated. the data is also attached, thank you.

clc
clear


JA = readmatrix('Acc.txt');
JA = JA-0.0154134; %takes it to zero
JA = detrend(JA);
tA = readmatrix('time.txt');
LVDT = readmatrix('LDVT.txt');
tL = readmatrix('t.txt');
Kistler = readmatrix('AccKis.txt');
Kistler = Kistler - 0.992151425253108; %takes it to zero
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, 4, 51);
Kistlersgf = sgolayfilt(Kistler, 4, 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')
xline(20, '-r', 'Filter Passband Frequency')
hold off;


[JAlpf, DF] = lowpass(JAsgf, 10, Fs, 'ImpulseResponse','iir');      % 20Hz Passband (Elliptic Filter)
[Kistlerlpf, DFK] = lowpass(Kistlersgf, 10, Fs, 'ImpulseResponse','iir');
% DF
figure
hold on;
plot(tA, JAlpf)
plot(tA, Kistlerlpf)
grid
xlabel('Time')
ylabel('Amplitude')
title('JA SG & Lowpass Filtered')
hold off

JAlpf = JAlpf*10^4;                       %to scale
Kistlerlpf = Kistlerlpf*10^4;
JV = cumtrapz(tA, JAlpf);                                           % Integrate To Get Velocity
JD = cumtrapz(tA, JV);                                              % Integrate To Get Displacement
KistlerV = cumtrapz(tA, Kistlerlpf);                                % Integrate To Get Velocity
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;

[Star Strider]

Just to be clear, what appears to be the problem? (I will likely get the same result you’re getting when I run your code, so I need to know what the problem is, since the code you posted appears to be quite similar to the earlier code I wrote.)

If the new acceleration data are similar to the previous acceleration data, using detrend will likely not change anything simply because there is no trend in the acceleration data. Subtracting the mean of the acceleration data vector from the acceleration vector would be helpful because it would remove a constant that would otherwise be integrated to produce a linear trend (that I do not see here) in the integrated data.

If the previous data is a guide, the acclerometer data (even when integrated) will not reproduce the LVDT data because the LVDT has a significantly different frequency response and a significantly lower sampling frequency than the accelerometer.

I will look at this in detail when I understand what I need to look for.

[Emily Keys]

Hi sorry, It was essentially the same issue where the acceleration data when double integrated to get displacement (as seen in the graph as the orange and purple lines) does not have a straight line before and after the dip that is the change in displacement. The pre and post shoulder (the time before and after the displacement dip can be seen) should essentially be 0 as there was no movement (as seen in the green line, the 'accurate' LVDT data which is used comparatively). It was a different experiment where the max displcament was 3mm instead of 24mm but the data seems to be a lot more skewed than it was for the data it was for the previous with the displacements having a lot of error. No worries if theres nothing that can be done and that is as accurate it can get. Thanks again.

[Star Strider]

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.

.