**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Plot spectra and spectrogram of acceleration data

35 views (last 30 days)

Show older comments

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

Mathieu NOE
on 14 Dec 2021

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
on 14 Dec 2021

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
on 14 Dec 2021

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)

Louise Wilson
on 14 Dec 2021

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
on 15 Dec 2021

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
on 16 Dec 2021

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
on 21 Dec 2021

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
on 22 Dec 2021

Edited: Louise Wilson
on 22 Dec 2021

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
on 3 Jan 2022

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
on 10 Jan 2022

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
on 11 Jan 2022

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
on 18 Jan 2022

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
on 19 Jan 2022

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

Mathieu NOE
on 16 Feb 2022

hello Louise and welcome back

I guess with datetick it should work :

all the best

### Accepted Answer

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
on 14 Dec 2021

Edited: Louise Wilson
on 14 Dec 2021

Bjorn Gustavsson
on 15 Dec 2021

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

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

### Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)