Filtering out noise or eigenfrequency in acceleration data

4 views (last 30 days)
Hello everybody,
For my master thesis I recorded some acceleration data at 4 measuring devices (see file attached, eacht column corresponds to one measuring device). The samling frequency is 300, so time step 1/300 = 0.0033s.
If I plot this graph (t against acc) I get lots of noise and high peaks. I guess they are because of the measuring device noise and some eigenfrequency of the specimen and the shaking table. I know from fft that the eigenfrequency of the shaking table is 50 Hz and from the specimen it is 6.85 Hz.
I tried many different filters but I did not find the correct one to clean out this high frequency acceleration peaks. Can you help me?
So far I did:
clc; clear;
load panels_umbro100_acc
acc = acc_g.*9.81; % because sampling data in in unit [g]
% create time vector
tStep = 0.00333; % Length of each time step
N = length(acc)*tStep;
t = 0:tStep:N;
t = t'; % column vector
t(end) = [];
dt = mean(diff(t)); % Average dt
fs = 1/dt;
% detrending
acc = detrend(acc);
% filter
N = 2; % quadratic
absolutevalue = 0.08;
[B,A] = butter(N,absolutevalue,'low');
acc2 = filter(B,A,acc); % filtered acceleration [m/s2]
% and I tried this filter to filter out fundamental frequency of specimen
% but what should I take as a q factor?
wo = 6.85/(300/2); % frequency to filter (6.85Hz) and sampling freq (300)
bw = wo/1; % wo / q factor (?)
[u,v] = iirnotch(wo,bw);
fvtool(u,v);
acc3 = filter(u,v,acc2); % filters acc2 by u and v
acc3_g = acc2./9.81; % filtered acceleration [g]
% plot for first measuring device
figure
figure1 = plot(t,acc3_g(:,1));
xlabel('t [s]')
ylabel('acc [g]')
grid on
The acceleration data of the first measuring device was on the shaking table so it should be the same as the graph attached.
I am really desperate and happy about any answer! Thanks a lot!

Accepted Answer

Mathieu NOE
Mathieu NOE on 9 Dec 2021
hello Lisa
do not be desperate ! we are here to help
I used your data in my favourite code for noise and vibration analysis. It will plot time , averaged fft spectra and spectrograms
The 50 Hz tone is a classicla electric perturbation from the mains. can be easily removed with a notch filter as you did (only my code does it with my own version of notch but this is not an issue).
then when you say you are in trouble with high frequency peaks , we can always use a low pass filters to remove them , but the question is how you know what belongs to the "true" signal and what is noise or perturbation ? the peaks can be valid data
the introduction of a low pass filter is not more complicated as the notch filter. We can use Butterworth filters (butter)
I added this also in my script
you can of course modify the filter type / order / cut off frequency according to your needs
all the best
code :
clc
clearvars
option_notch = 1; % 0 = without notch filter , 1 = with notch filter
fc_notch = 50; % notch freq
option_LPF = 1; % 0 = without low pass filter , 1 = with low pass filter
fc_lpf = 75; % LPF cut off freq
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
load panels_umbro100_acc
% detrending
acc_g = detrend(acc_g);
signal = acc_g.*9.81; % because sampling data in in unit [g]
[samples,channels] = size(signal);
% create time vector
Fs = 300;
dt = 1/Fs;
time = (0:samples-1)*dt;
%% notch filter section %%%%%%
% y(n)=G*[x(n)-2*cos(w0)*x(n-1)+x(n-2)]+[2*p cos(w0)*y(n-1)-p^2 y(n-2)]
% this difference equation can be converted to IIR filter numerator /
% denominator
if option_notch ~= 0
w0 = 2*pi*fc_notch/Fs;
p = 0.995;
% digital notch (IIR)
num1z=[1 -2*cos(w0) 1];
den1z=[1 -2*p*cos(w0) p^2];
% now let's filter the signal
signal = filtfilt(num1z,den1z,signal);
end
%% low pass filter section %%%%%%
if option_LPF ~= 0
w0_lpf = 2*fc_lpf/Fs;
% digital notch (IIR)
[b,a] = butter(2,w0_lpf);
% now let's filter the signal
signal = filtfilt(b,a,signal);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512; %
OVERLAP = 0.95;
% spectrogram dB scale
spectrogram_dB_scale = 80; % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums
% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
for ck = 1:channels
newsignal(:,ck) = decimate(signal(:,ck),decim);
Fs = Fs/decim;
end
signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;
%%%%%% legend structure %%%%%%%%
for ck = 1:channels
leg_str{ck} = ['Channel ' num2str(ck) ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),plot(time,signal);grid on
title(['Time plot / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(leg_str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);
% apply A weigthing if needed
if option_w == 1
pondA_dB = pondA_function(freq);
sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
my_ylabel = ('Amplitude (dB (A))');
else
my_ylabel = ('Amplitude (dB (L))');
end
figure(2),plot(freq,sensor_spectrum_dB);grid on
df = freq(2)-freq(1); % frequency resolution
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
legend(leg_str);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
[sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg)); % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% apply A weigthing if needed
if option_w == 1
pondA_dB = pondA_function(fsg);
sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
my_title = ('Spectrogram (dB (A))');
else
my_title = ('Spectrogram (dB (L))');
end
% saturation of the dB range :
% saturation_dB = 60; % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2+ck);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid on
df = fsg(2)-fsg(1); % freq resolution
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
xlabel('Time (s)');ylabel('Frequency (Hz)');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function pondA_dB = pondA_function(f)
% dB (A) weighting curve
n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
pondA = n./r;
pondA_dB = 20*log10(pondA(:));
end
function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).
% Linear averaging
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,channels);
s_tmp((1:samples),:) = signal;
signal = s_tmp;
samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
% compute fft with overlap
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
% % for info is equivalent to :
% noverlap = Overlap*nfft;
% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows
% main loop
fft_spectrum = 0;
for i=1:spectnum
start = (i-1)*offset;
sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only
end
fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end
  2 Comments
Mathieu NOE
Mathieu NOE on 9 Dec 2021
hello Lisa
no problem - my pleasure
do not hesitate to ask if needed
all the best
M

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by