Plot spectra and spectrogram of acceleration data

Question

Louise Wilson 2021-12-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1610145-plot-spectra-and-spectrogram-of-acceleration-data

评论： Mathieu NOE 2022-2-16

example.mat

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 个评论
显示 13更早的评论隐藏 13更早的评论

Mathieu NOE 2021-12-15

hello Louise

seems the wav file is not appearing right now (still the original mat file) - can you try again ?

I don't know excatly what kind of recorder you are using but if the y "true" scale matters in your spectral analysis, there are a few points to check

did I use the sensor's sensivity in my recorder set up ? I had a look at your sensor datasheet :

Output sensitivity:

Accelerometer : 10 V/g

Hydrophone -162 dB re 1.0 V/uPa

in my recorder file format , does the data appears in volts or directly converted into the engineering units (according to sensivity given above - do not forget to pay attention if you use a separate sensor amplifier with a non unity gain)
pay attention to what happens when you export your data in wav format. the wav format option will probably apply a scaling factor to match the +/- 1 max amplitude format of wav files (beyond these limits the signal will be clipped => distorted). So if you have a signal which is in a larger range (like +/- 10 range) it will be rescaled down to match +/- 1 max range. If you are not aware of this process , you might think when analysing the wav file that the signal has lost a factor 10 of amplitude. If you have a signal that is already within the limts of +/-1 amplitude range , the export in wav format should not modify the dynamic range, but it's alaways good practive to double check. Sometimes the wav export process will tell you by how much the original signal amplitude has been rescaled to avoid clipping, but that's specific to recorder's / analyser's vendors.

Mathieu NOE 2021-12-21

在 MATLAB Online 中打开

hello Louise

I tooke the initial portion of your code to apply the sensivity factors

the rest of the code is my usual "averaging fft" and spectrogram. i believe there is only a vertical shift (in dB) between your spectrum amplitude and mine but the spectral content is the same (major peak around 4kHz)

in y plot below o dB = 1 m/s² engineering unit

the spectrograms are similar

clc
clearvars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 
Fs=48000;
signal=example.xacc{1}; %read data
dt = 1/Fs;
[samples,channels] = size(signal);
time =  (0:samples-1)*dt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 512;    % 
OVERLAP = 0.75;
% 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

Louise Wilson 2021-12-22

编辑：Louise Wilson 2021-12-22

在 MATLAB Online 中打开

Thanks again for your dedication to helping me solve this problem Mathieu :)

Before I heard back from you I had tried again a different way, and got the following...

% Input variables
fs = 48000;
nfft=1624;
noverlap=round(0.5*nfft);
w=hanning(nfft);
% Spectra
subplot(2,1,1)
[sensor_spectrum, freq] = pwelch(example.xacc{1,1},w,noverlap,nfft,fs);
sensor_spectrum_dB = 10*log10(sensor_spectrum);
plot(freq,sensor_spectrum_dB);grid on; %or semilogx()
ylabel('Power Spectral Density (dB re (1ms^-^2)^2 / Hz)');
xlabel('Frequency');
% Spectrogram
subplot(2,1,2)
[S,F,T,P] = spectrogram(example.xacc{1,1},w,noverlap,nfft,fs);
surf(T,F,10*log10(P),'edgecolor','none');
colormap(jet);                                                             
axis tight;
view(0,90);
colorbar
ylabel('Frequency (Hz)','fontsize',12, 'FontName', 'Arial');
ylabel(colorbar,'Power Spectral Density (dB re (1ms^-^2)^2 / Hz)','Rotation',270,'VerticalAlignment','bottom','FontSize',12);
set(gca,'TickDir','out'); 

As before, there is a shift in the dB values between my and your outputs. So I guess there is something wrong with my code and I will spend some time interpreting yours. But, do you agree the resolution of the spectrogram looks better here? I had thought this was down to overlap and nfft size but with the same parameters I can't get a similar output with your code?

Mathieu NOE 2022-1-3

在 MATLAB Online 中打开

hello Louise

happy new year for you and your familly !!

I will have again a look at the dB shift - probably tomorrow

for the spectrogram look , the frequency resolution df only depends on Fs and nfft (df = Fs/nfft)

so if both codes do share the same Fs and nfft values there should be no difference in frequency resolution; you can add a second line in your code to compute it :

[sensor_spectrum, freq] = pwelch(example.xacc{1,1},w,noverlap,nfft,fs);
df = freq(2)-freq(1); % freq resolution 

the time resolution can be improved by increasing the overlap , a higher overlap means we create more intermediate spectrum slices in the spectrogram so the visual aspect seems better; yiu can play with the overlap factor to see that effect; If you need a refined time resolution you can go up to 95% of overlap (50 or 75% is my usual start value);

also the other parameters that may create a different visual rendering of the spectrogram is the colormap (here we have both jet), and aslo in my code I can select the dB range so that what is below the lower limit (= upper limit minus the selected dB range) is discarded; this way I can really zoom into a given dB range if I need to; can be interesting if you have two frequencies of similar amplitudes but you need to make sure which one is higher in amplitude vs the other one;

all the best

Mathieu NOE 2022-1-11

在 MATLAB Online 中打开

hello Louise

how are you ?

I noticed that maybe there was something wrong when you apply the sensivity factors.

for me, the 9.81 factor should be applied at the numerator side and not in the denominator

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

(you can also see that if I convert the sensivity of the accel from 10V/g to 10/9.81 ms-2 / V , so at the end you would divideyour data in volts by 10/9.81 , which is the same as multiply by 9.81 and divide by 10)

otherwise , in your code , you could apply the same trick I used in my own version to force the data to lie within a given dB range

this should bring you the same spectrogram plot as in my code

otherwise maybe there are a ew differences between how spectrogram amplitudes are computed between the older specgram function I am still using and the newer spectrogram version. haven't really taken the time to compare both functions yet to be honest.

your code a bit updated :

% Input variables
fs = 48000;
nfft=512;
noverlap=round(0.75*nfft);
w=hanning(nfft);
% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
% Spectra
figure
subplot(2,1,1)
[sensor_spectrum, freq] = pwelch(example.xacc{1,1},w,noverlap,nfft,fs);
sensor_spectrum_dB = 10*log10(sensor_spectrum);
plot(freq,sensor_spectrum_dB);grid on; %or semilogx()
ylabel('Power Spectral Density (dB re (1ms^-^2)^2 / Hz)');
xlabel('Frequency');
% Spectrogram
subplot(2,1,2)
[S,F,T,P] = spectrogram(example.xacc{1,1},w,noverlap,nfft,fs);
    % update MN
    sg_dB = 10*log10(P);
    % saturation of the dB range : 
    min_disp_dB = round(max(sg_dB,[],'all')) - spectrogram_dB_scale;
    sg_dB(sg_dB<min_disp_dB) = min_disp_dB;
surf(T,F,sg_dB,'edgecolor','none');
colormap(jet);                                                             
axis tight;
view(0,90);
colorbar
ylabel('Frequency (Hz)','fontsize',12, 'FontName', 'Arial');
ylabel(colorbar,'Power Spectral Density (dB re (1ms^-^2)^2 / Hz)','Rotation',270,'VerticalAlignment','bottom','FontSize',12);
set(gca,'TickDir','out'); 

and the output (for the first wav file you supplied)

Louise Wilson 2022-2-14

编辑：Louise Wilson 2022-2-14

在 MATLAB Online 中打开

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?

Mathieu NOE 2022-2-16

hello Louise and welcome back

I guess with datetick it should work :

Date formatted tick labels - MATLAB datetick - MathWorks France

all the best

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bjorn Gustavsson 2021-12-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1610145-plot-spectra-and-spectrogram-of-acceleration-data#answer_854875

编辑：Bjorn Gustavsson 2021-12-14

在 MATLAB Online 中打开

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 个评论
显示 1更早的评论隐藏 1更早的评论

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 2021-12-21

在 MATLAB Online 中打开

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]);

请先登录，再进行评论。

Plot spectra and spectrogram of acceleration data

15 个评论
显示 13更早的评论隐藏 13更早的评论

采纳的回答

3 个评论
显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Plot spectra and spectrogram of acceleration data

15 个评论 显示 13更早的评论隐藏 13更早的评论

采纳的回答

3 个评论 显示 1更早的评论隐藏 1更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

15 个评论
显示 13更早的评论隐藏 13更早的评论

3 个评论
显示 1更早的评论隐藏 1更早的评论