Plot spectra and spectrogram of acceleration data

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

[Mathieu NOE]

hello Louise

welcome back !

the code example I provided in the other post (this question) does not make / need any assumption about where the measurement is done below or above sea level...

the only difference may arise from sensor's sensivity (which is not included in my code) , dB reference levels (they differ if you speak vibration or noise dB spectra) and sampling rate

if your spectra looks like having too much level (in dB ) have you checked if the time data has the correct scaling (engineering unit ? ) what is the sensor used and did you convert the raw (voltage) signal into a physical unit signal ?

if the frequency range seems too broad , I would ask myself if we are sure the signal rate is actually 48 kHz or if any in between processing has been done to modify the original sampling rate (downsampling ?) ..

hope this helps

all the best

[Mathieu NOE]

I didn't pay attention to this info and @Bjorn Gustavsson is right about that

if we assume Fs = 48 kHz (as in your code) , the signal duration is 5 seconds and not 5 minutes

if the later is true then Fs = 800 Hz as confirmed by Bjorn

this will rescale all your frequency values by a factor of 1/60 which may now be more in line with your expectations

for the y amplitude look for sensor senvity and amplifier's gain that are not included in your / my code yet.

[Mathieu NOE]

BTW, I am also intriguished by the fact that your are analysing boat sound with accelerometers - that's a bit unusual... would have thougt you are using hydrophones... unless this is a "home made" hydrophone !

also a sampling of 800 Hz seems to me too low now ... as boat noise covers a much broader range than that (cavitation of propellers goes in the x kHz range, machinery noise mostly below 1 kHz)

Revolutionary Propeller Technology Reduces Underwater Radiated Noise (maritime-executive.com)

[Louise Wilson]

Thanks Mathieu!

I have now updated the question so you can access the raw data and see how I read and calibrated the data. The sampling rate of 48 kHz is reported when I read in the data and this is also the sampling rate I programmed the recorder with.

The sensor used is a low frequency vector sensor made by Wilcoxon (VS-301). The flat frequency response is 3 Hz to 2 kHz.

I think that the reference level also varies underwater?

Yes, I have used hydrophones previously, but the reason to work with an accelerometer is because we are interested in the effects of the boat sound on invertebrates and fish, who hear by particle motion :)

[Mathieu NOE]

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.

[Louise Wilson]

Hi Mathieu, can you see the file here? https://drive.google.com/file/d/1vC8bE4NbvG8iVFwrqBr2LccW3c7N9G7w/view?usp=sharing

The recorder is a SoundTraps recorder, the sensitivity and attenuation for this is the bd_sen and bd_att. Then those for the accelerometer/vector sensor are below. I am not using the hydrophone in this instance. The recorder outputs a .wav file every 5 minutes which has four channels (X,Y,Z,hydrophone).

I believe that in the file, the units are in Volts, but then I calibrate them and detrend in the code by applying the above sensitivities and attenuation? Now the units are in m/m/s.

[Mathieu NOE]

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]

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]

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

[Louise Wilson]

Hi Mathieu,

Happy New Year to you and yours also :)

Thanks a lot for your perseverence with this question!

I compared the frequency resolution of our approaches and for both - the value is 29.5567. So that is reassuring! I see now that this is the value you report in the title of the spectra plot :)

Using 75% overlap and a colorbar limit of -120 to -50 we get quite different outputs for the spectrogram so I guess this is down to the dB shift?

In your plot, you can see the approach of a second boat pass at the end of the file (> ~250 seconds) (which I know to be the case) whereas this is not evident in mine, showing some detail has been missed.

Mathieu:

inexperienced Louise:

[Mathieu NOE]

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]

Thanks so much Mathieu!

I can't explain it but I originally calculated the sensitivity the way you have corrected. But, I thought that the way I have done it is correct as I used two methods to calculate the acceleration. I used an acceloremeter (data we are working with here) and also a hydrophone array where I estimate the acceleration using other calculations. The values only make sense (are similar) if I calculate the sensitivity as in my original script here. And in this case the values are close to what I would expect from other published data.

Here in orange is the data from the accelerometer (more sensitive) compared with the array (blue) which makes sense as the array method is an approximation and not likely to be as accurate. Using the different approach for sensitivity we get values in the range 0 - 0.5 m/s/s which doesn't really make sense, it is too fast?

[Mathieu NOE]

hello Louise

i don't want to add any confusion here ...there must be for sure good reasons why your method works better than what I suggested

from my own experience of vibration measurement, I (still) stick to the fact that if I have a 10 V / g accel , I have to divide by signal in volts by 10 if I want the result in g , and then multiply by 9.81 if I want the result in m/s² , which is basically the same as divide volts by a new sensivity S2 = 10/9.81 (rounded to 1 V / m/s²)

all the best

[Louise Wilson]

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]

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]

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

[Louise Wilson]

Thanks a lot Bjorn. I am still not sure if the sampling rate is wrong? I programmed the recording device to sample at a rate of 48 kHz, but maybe there is something here I don't understand. I can see that your figure makes more sense in terms of S. However, the frequency range of the device is 50 - 2000 Hz so I would have expected it to be broader than 0-400 Hz? I have now updated the question to included the .wav file and the processing I applied to calibrate the data for the recording devices used and convert the units.

[Bjorn Gustavsson]

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]

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