Plot spectra and spectrogram of acceleration data

35 views (last 30 days)
I have acceleration data for three planes (X,Y,Z) which is recorded using an accelerometer connected to a recorder (referred to as board).
An example .wav file is here: https://drive.google.com/file/d/1vC8bE4NbvG8iVFwrqBr2LccW3c7N9G7w/view?usp=sharing
I would like to plot the spectra and spectrogram of each plane.
I tried using this question as a starting point to work on the outputs for one plane. However, I assume this question is specific to terrestrial measurements of acceleration and my data was collected using an accelerometer deployed in the sea at a depth of ~ 8m. The example file is 5 minutes long and contains the sound of a boat passing by.
filename='805367873.210811101406.wav'; %wav file (5 minutes long)
%Board and vector sensor sensitivities
bd_sen = 1.363; % 5.7-3 dB (1.363 is Volts/unit)
bd_att = 4 % 12.04 dB
vs_sen_x = 10; %Sensitivies for vector sensor #166 (V/g)
vs_sen_z = 10.2;
[D, fs] = audioread(filename);
%fs is 48000
%Calibration values
x_sen = (bd_att*bd_sen/(9.81*vs_sen_x)); %ms^-2/unit
z_sen = (bd_att*bd_sen/(9.81*vs_sen_z)); %ms^-2/unit
%Apply calibration and detrend
example.xacc{1} = detrend(D(:, 1)*x_sen); %x
example.zacc{1} = detrend(D(:, 3)*z_sen); %z
%FFT parameters
nfft=1624;
noverlap=round(0.5*nfft);
w=hanning(nfft);
spectrogram_dB_scale=80;
samples=length(signal);
Fs=48000;
% Averaged FFT spectrum
[sensor_spectrum,freq]=pwelch(signal,w,noverlap,nfft,fs);
sensor_spectrum_dB=20*log10(sensor_spectrum) %maybe this should be 10*log10(spectrum)?
figure(1),semilogx(freq,sensor_spectrum_dB);grid
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');
%The values here seem too large, and the frequency range is broader than
%I would have expected. I would have expected an upper limit of around 3000
%Hz.
%Spectrogram
[sg,fsg,tsg] = specgram(signal,nfft,Fs,w,noverlap);
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
xlabel('Time (s)');ylabel('Frequency (Hz)');
  15 Comments

Sign in to comment.

Accepted Answer

Bjorn Gustavsson
Bjorn Gustavsson on 14 Dec 2021
Edited: Bjorn Gustavsson on 14 Dec 2021
Provided that your data actually spans 5 minutes your sampling-frequency is not 48 kHz, but rather 800 Hz. That changes the output a bit:
fs = (numel(example.x)-1)/5/60;
nfft=1624;
noverlap=round(0.5*nfft);
w=hanning(nfft);
[S,F,T,P] = spectrogram(example.x,w,noverlap,nfft,fs);
pcolor(T,F,log10(abs(S))),shading flat,colormap(turbo),colorbar
if you have band-passed downconverted data this should be adapted, but for that you'll have to explain what is done to the data.
HTH
  3 Comments
Louise Wilson
Louise Wilson on 21 Dec 2021
Thanks Bjorn. I didn't realise that I made a mistake uploading the data, hence the fs confusion! So, this question was doomed from the start... Once I resolved that I figured this out with help from your code, although I changed things slightly:
fs = 48000;
nfft=1624;
noverlap=round(0.5*nfft);
w=hanning(nfft);
[S,F,T,P] = spectrogram(signal,w,noverlap,nfft,fs);
surf(T,F,10*log10(P),'edgecolor','none');
colormap(jet);
axis tight;
view(0,90);
colorbar
caxis([-200 -65]);

Sign in to comment.

More Answers (0)

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!