Plot spectra and spectrogram of acceleration data

40 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
Louise Wilson
Louise Wilson on 14 Feb 2022
Edited: Louise Wilson on 14 Feb 2022
Hi Mathieu,
I have been using your code to plot some acceleration data. The data is in 5 minute chunks (because that's how the recorder works). In our examples so far we just worked on one 5 minute file.
So, to plot all of the data together (many five minutes files), I concatenated all of the data to produce a 102000425x1 array, and plot this as below:
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg)); % NB : X=fft(x.*hanning(N))*4/N; % hanning only
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');
%time, frequency, colour
axis('xy');%colorbar('vert');
colorbar
colormap jet
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)');
This looks good, but I'd like to have a more informative x axis. So I considered that I could add the start time of the first recording to create an array of datetimes:
tsg_dt(1)=datetime(2021,11,8,9,51,0);
for i=2:height(tsg)
tsg_dt(i)=tsg_dt(1)+seconds(tsg(i));
end
But, this doesn't work with imagesc or surf.
imagesc(tsg_dt,fsg,sg_dBpeak);colormap('jet');
%Error using image
%Image XData and YData must be vectors.
surf(tsg_dt,fsg,sg_dBpeak,'FaceColor','interp',...
'EdgeColor','none','FaceLighting','phong');
view(0,90)
%"works" but does not look right after changing the view (maybe this is the
%main problem...?)
Do you have an idea for how else I could do this?

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
Bjorn Gustavsson
Bjorn Gustavsson on 15 Dec 2021
The question for you to resolve first is how many samples do you have from a time-period that is how long and does that match up with your sampling-frequency-setting? If you have base-band samples at 48 kHz then you should have 48000 samples per second (obviously). Does that match up with the length of your data-set, 5 minutes, or 5 seconds? Are fs and Fs identical?
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!