How to calculate PSD estimate of acoustic data in 1Hz bins? pwelch?

Hi everyone,

I have several files of 118.5 second duration that I would like to calculate the PSD of. I want to do this in 1Hz frequency bins, with Hanning window and 50% overlap, and have considered using pwelch.

First I calibrate the data and then have tried to work out how to use pwelch to no avail so far:

path=('C:\Users\datafolder');

d=dir(fullfile(path, '*.wav')); %list all .wav files in path folder

files=length(d); %number of files in folder

for i=1:files    %for each file
    disp(d(i).name);
    filename=fullfile(path, d(i).name); %get full filename
    [xbit, fs]=audioread(filename); %and then sampled data and sampling 
                                    %freq (fs) from audioread function
    cal=-176.2; 
    cal=power(10, cal/20); %calibration correction factor
    calxbit=xbit*cal; %apply calibration
end


[pxx,f]=pwelch(calxbit, window, 0.5, [], 144e3);

How do I specify that I want to use 1Hz bins? My data is sampled at 144kHz.

Thanks!!

[Michael Ladegaard]

One of your problems stems from the calibration value, which I assume is in the unit dB re V/µPa. First, you simply need to change the sign of this calibration value, so that it becomes 176.2 dB re µPa/V. If the signal 'xbit' is in Volts, then you are left with units of pressure after you multiply the 'calibration correction factor' to your signal:

    cal=-176.2; % dB re V/µPa
    cal=-cal; % dB re µPa/V    
    cal=power(10, cal/20); %calibration correction factor, µPa/V
    calxbit=xbit*cal; %apply calibration, Unit is µPa.

Secondly, it looks like you want 50% overlap, but the overlap input should not be a ratio such as 0.5, but an integer such as round( 0.5*length(window) ).

Lastly, if you want to compute the PSD for each wave file, you need to move the pwelch-line into the for loop (see code below).

If you instead want to perform the pwelch routine for all the audio data, you need to concatenate the 'calxbit'-signals - currently you just overwrite the 'calxbit' and 'pxx' variables in every loop.

path=('C:\Users\datafolder');
d=dir(fullfile(path, '*.wav')); %list all .wav files in path folder
files=length(d); %number of files in folder

% pwelch parameters
nfft = 1024 ; % FFT size
window = hann(nfft); % define the window
noverlap = round(nfft/2) ; % an integer

for i=1:files    %for each file
    disp(d(i).name);
    filename=fullfile(path, d(i).name); %get full filename
    [xbit, fs]=audioread(filename); %and then sampled data and sampling 
                                    %freq (fs) from audioread function
    xbit = detrend(xbit) ; % remove DC
    cal=-176.2; % dB re V/µPa
    cal=-cal; % dB re µPa/V    
    cal=power(10, cal/20); %calibration correction factor. µPa/V
    calxbit=xbit*cal; %apply calibration, Unit is µPa if xbit is in units of Volts.
    
    [pxx,f]=pwelch(calxbit, window, noverlap, nfft, fs,'psd'); % 
    PXX{i} = pxx ; % Save each estimate in a cell so it is not overwritten in each loop
    F{i} = f ;
end

[Mathieu NOE]

regarding the fft frequency resolution (df) question : How do I specify that I want to use 1Hz bins? My data is sampled at 144kHz

remember that df = fs/nfft

so in order to get a 1 Hz resolution it requires nfft = fs. Seems to me quite huge as we know fs is 144kHz. It may take some time to compute such a long fft....

also , a minor improvement can be to apply the calibration factor on Pxx instead of making very long array multiplications when you apply the correction on the time signal (if the time signal does not need conversion in microP , just the Pxx for display)

[Louise Wilson]

Thanks @Mathieu NOE - that is indeed what I did in the end. It takes a long time but is the resolution I need. I also applied the calibration factor on pxx, by converting pxx to dB. :)

[Louise Wilson]

Thanks @Michael Ladegaard - sorry, I did not realise your comment was also recent :). In the end I did account for the three points you mention, but I am sure your solution is more elegant! I will refer back to it. Thanks a lot :)