Azzera filtri
Azzera filtri

Frequency based on acceleration data

22 visualizzazioni (ultimi 30 giorni)
MARLIANA HAFIZA
MARLIANA HAFIZA il 14 Apr 2023
Commentato: Mathieu NOE il 12 Mag 2023
Hi,
I want to find the frequency during testing based on the acceleration data collected using accelerometer. I already try several method using MATLAB but I'm not sure if mine is correct. Different method give me different frequency value. For the info, the machine frequency is 6.67 Hz. I attach the code and file for reference.
clc;
clear all;
close all;
Time = xlsread('Book5','A1:A585');
Acceleration = xlsread('Book5','B1:B585');
n = length(Acceleration);
Ts = Time(2)-Time(1);
Fs = 1/Ts;
NFFT = 2^nextpow2(n);
Y = fft(Acceleration,NFFT)/n;
f = Fs/2*linspace(0,1,NFFT/2+1);
Iv = 1:length(f);
figure(1)
plot(f,abs(Y(Iv)))
xlabel('Frequency (Hz)')
ylabel('Amplitude (m)')
num=csvread("Book5.csv");
Fs = 1/mean(diff(num(:,1)));
Fn = Fs/2;
Tr = linspace(num(1,1), num(1,end), size(num,1));
Dr = resample(num(:,4), Tr);
Dr_mc = Dr - mean(Dr,1);
FDr_mc = fft(Dr_mc, [], 1);
Fv = linspace(0, 1, fix(size(FDr_mc,1)/2)+1)*Fn;
Iv = 1:numel(Fv);
figure(13)
plot(Fv, abs(FDr_mc(Iv,:))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Amplitude')
Thank you.

Risposte (1)

Mathieu NOE
Mathieu NOE il 14 Apr 2023
hello
I tried to do my best finding the fft length and overlap that gives me the best info - as your signal is quite unstationnary
this is what I could tell from your data
Spectrogram : look at the red area as this is the dominant frequencies (depend of time)
the amplitude (color) is expressed in dB (log scale of amplitude)
Averaged FFT spectrum , basicaly the time average of the spectrogram , here you lose the information of time dependance but you can probably find your frequencies (no 6.67 Hz found ! )
Code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename = ('Book5.xlsx');
data = readmatrix(filename);
time = data(:,1);
signal = data(:,2);
dt = mean(diff(time));
Fs = 1/dt; % sampling rate
[samples,channels] = size(signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 200; %
OVERLAP = 0.95; % must be between 0 and 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
  3 Commenti
Mathieu NOE
Mathieu NOE il 17 Apr 2023
in terms of signal processing, you can decide wheter you make
  • one fft computation with the fft length equal to the signal length (and in this case there is of course no averaging
  • or you split your signal into smaller chunks (buffers) and you do the fft on each buffer and then you average the spectra ; this reduce the impact of spurious / non correlated noises on your resulting spectrum. You can increase the number of averages by introducing some overlap between buffers. For spectrograms, having / increasing the overlap makes the rendering more progressive (smooth) at the expense of more fft computations
  • also the fft spectrum frequency resolution df depends of the fft length nfft and sampling rate (Fs) :
  • df = Fs / nfft, so for a given Fs, if you want more averages you will have to reduce nfft (and maybe use a higher overlap), but globally you cannot have a the same time a very small frequency resolution and lot of averages. It's all about finding the best compromise depending of your signal characteristics (stationnary or not)
Mathieu NOE
Mathieu NOE il 12 Mag 2023
hello again
did my submission fullfil your request ?

Accedi per commentare.

Tag

Community Treasure Hunt

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

Start Hunting!

Translated by