Matlab Spectrogram: How to convert the dB values to SPL?

Guys,

So I have to plot this spectrogram, but of course when I did the magnitude scale came out in -dB. Is there some way to change this to SPL, Pa?

The following is the code used;

if true
   [x, Fs] = wavread ('C:\Users\User\Desktop\3020 Recordings_Final\C#\Straight Legs\Jawari Design 1 (Staright Legs)\Apple Ma\No Fret_d1am');%Reads file into Matlab
D = x(:,1);                                     % Extracting channel 1 data
% plotting of the spectrogram
figure (21)
spectrogram(D, 1024, 3/4*1024, [], Fs, 'yaxis')
h = colorbar;
set(h, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(h, 'Magnitude (dB');
xlabel('Time, s', 'fontsize', 16)
ylabel('Frequency, Hz', 'fontsize', 16)
title('Spectrogram of the signal', 'fontsize', 20)

end

Also I have the calibration factor from using a pistonphone; 85.91Pa

My problem is, I don't know what Matlab is doing within that function when it calculates the dB value.

Any help will be much appreciated, thanks!

[Joaquin]

Hello,

I have exactly the same question. If you have managed to solve it, I would appreciate it if you could tell me how, because I can't get the colorbar in the spectrogram to reflect SPL in Pa or dB.

Thank you so much,

Joaquin

Answer 1

[Mathieu NOE]

hello

I admit it's not always evident to know how n normalize spectrogram data , especially when now in the new matlab releases the main code is hidden.

I ended up doing my own spectrogram code where I know what I do . I was just a bit lazzy so I took only one window (hanning) by default and it's corresponding amplitude correction coefficient

the code belows simply compute the spectrogram data and you can correct for the calibration factor , assuming you have first read the calibration wav file (pistonphone), then once the scale _factor is known (see to the code ) you must hard code that value for the next wav files

here for example a created a sine wave with 0.23 amplitude (this could be what you have in your wav calibration file) , that is supposed to correspond to the 85.91Pa level

therefore here the scale factor is obtained in this line : scale_factor = 85.91./max(abs(D));

full code :

%%%%%%%%%%%%%%%%%%%%%%%
% dummy data
Fs = 2e3;
dt = 1/Fs;
samples = 10000;
t = (0:samples-1)'*dt;
x = 0.23*sin(2*pi*100*t); % sine wave of amplitude whatever , obtained by calibration (pistonphone)
%%%%%%%%%%%%%%%%%%%%%%%

D = x(:,1);                                     % Extracting channel 1 data

% scale factor : max amplitude read from wav file correspond to : pistonphone; 85.91Pa
scale_factor = 85.91./max(abs(D));

% plotting of the spectrogram
nfft = 1024;
overlap = 3/4;
[S,F,T] = myspecgram(D, Fs, nfft, overlap); 

% apply scale factor 
S = S*scale_factor;

figure (21)
imagesc(T,F,S)
set(gca,'YDir','normal')
h = colorbar;
set(h, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(h, 'Magnitude');
xlabel('Time, s', 'fontsize', 16)
ylabel('Frequency, Hz', 'fontsize', 16)
title('Spectrogram of the signal', 'fontsize', 20)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function  [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal  (example sinus amplitude 1   = 0 dB after fft).
%   signal - input signal, 
%   Fs - Sampling frequency (Hz).
%   nfft - FFT window size
%   Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,1);
    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_specgram = [];
    for ci=1:spectnum
        start = (ci-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        fft_specgram = [fft_specgram abs(fft(sw))*4/nfft];     % X=fft(x.*hanning(N))*4/N; % hanning only 
    end

% 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_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;

% time vector 
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(nfft/2))/Fs;

end