Filtering out noise or eigenfrequency in acceleration data

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!

Answer 1

[Mathieu NOE]

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

[Lisa Oswald]

Dear Mathieu NOE!

Wow thank you very much for your time, effort and help! I need a bit to understand everything what you did. Will do that the next days.

Thanks a lot!

[Mathieu NOE]

hello Lisa

no problem - my pleasure

do not hesitate to ask if needed

all the best

M