Define main frequencies of accelerometer signals

Dear all,

I have some acceleration data, given by ambient vibrations, as the ones in the attached .txt file (X and Y on the same accelerometer), with a frequency sampling of 200. The main problem is due to a noise that does not allow to recognize the main frequencies of the signal and that pushes to perform some filtering actions (e.g., pass band, low or high pass). Anyway, I'm trying to understand what are the main frequencies of the signals and how to select the values to filter, also considering that PSD does not return clear peaks (an attempt was made with a low pass filter usin 20 Hz as frequency, for which I attached 2 figures, X and Y).

Could anyone help me?

Thanks in advance

[Image Analyst]

How about some screenshots? PNG files instead of fig or use the green frame icon to insert into the body of your message.

[Sergio RUGGIERI]

Ok, see the attached PNG files

Answer 1

[Mathieu NOE]

hello Sergio

if you are interested in the low frequency range , my liitle code can help you

the quantization noise is quite high but we still can see some emergence of 3 peaks around 14 , 17 and 20 Hz

try this

clc
clearvars

option_notch = 0; % 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 = 50;       % LPF cut off freq
    N_lpf = 2;           % filter order

option_HPF = 0; % 0 = without high pass filter  , 1 = with high pass filter
    fc_hpf = 150;       % HPF cut off freq
    N_hpf = 2;           % filter order

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

signal = importdata('dataAcc.txt');

[samples,channels] = size(signal);

% create time vector
Fs = 200;
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(N_lpf,w0_lpf);
    % now let's filter the signal 
    signal = filtfilt(b,a,signal);
end

%% high pass filter section %%%%%%
if option_HPF ~= 0
    w0_hpf = 2*fc_hpf/Fs;
    % digital notch (IIR)
    [b,a] = butter(N_hpf,w0_hpf,'high');
    % now let's filter the signal 
    signal = filtfilt(b,a,signal);
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512;    % 
OVERLAP = 0.75;

% dB scale
dB_scale = 40;  % 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);
mm = 5*ceil(max(sensor_spectrum_dB,[],'all')/5); % 5 dB rounded max level
ylim([mm-dB_scale mm]);
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))) - 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

sound(signal,Fs);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

[Sergio RUGGIERI]

Dear Mathieu,

thanks for your helpful code. It works fine.

Kind Regards.

Sergio

[Mathieu NOE]

hello Sergio

my pleasure !