Plot spectra and spectrogram of acceleration data

39 次查看(过去 30 天)
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 个评论
Louise Wilson
Louise Wilson 2022-2-14
编辑:Louise Wilson 2022-2-14
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?

请先登录,再进行评论。

采纳的回答

Bjorn Gustavsson
Bjorn Gustavsson 2021-12-14
编辑:Bjorn Gustavsson 2021-12-14
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 个评论
Bjorn Gustavsson
Bjorn Gustavsson 2021-12-15
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 2021-12-21
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]);

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Spectral Measurements 的更多信息

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by