Clear Filters
Clear Filters

How to use pwelch for PSD?

15 views (last 30 days)
Louise Wilson
Louise Wilson on 26 Jul 2019
Answered: emehmetcik on 8 Nov 2019
I would like to use pwelch to perform PSD on my long term acoustic data. The data is 120 seconds long, sampled at 144kHz.
How do I decide what to set nfft as? Is it correct that I have set it as the length of my sampled data?
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; %for ST5101, adapt this to use switch and have diff cal for
%diff ST used, should it be negative or not?
cal=power(10, cal/20); %calibration correction factor
calxbit=xbit*cal; %apply calibration
end
%% pwelch
%[Pxx,f]=pwelch(x, window, noverlap, nfft,fs)
%'x' is signal
%'window' divides the signal into segments.
%Default window length for PamGuide is N=Fs, yielding a frequency
%resolution of 1Hz.
%'noverlap' must be a positive integer less than length of window. If you do
%not specify noverlap the default number of overlapped samples is 50% of
%the window length.
%'nfft' number of discrete Fourier transform points to use in the DFT
%estimate.
nfft=length(calxbit);
%In our case, x is a vector (for each file).
[pxx,f]=pwelch(calxbit, segmentLength, 0.5, nfft, fs);
%[] is nfft, how to decide what to put here?
plot(f, 10*log10(pxx));
When I run this code, I get the error:
Warning: Integer operands are required for colon operator when used as index.
> In welch>localComputeSpectra (line 316)
In welch (line 99)
In pwelch (line 162)
...multiple times!
Thanks for your help!

Answers (1)

emehmetcik
emehmetcik on 8 Nov 2019
The way you call "pwelch" function is erroneous:
[pxx,f]=pwelch(calxbit, segmentLength, 0.5, nfft, fs);
When "pwelch" is called in this manner it expects the inputs to be as follows;
[pxx,w] = pwelch(x, window, noverlap, nfft, fs)
"noverlap" must be an integer less than nfft.
See the help page for pwelch for details.
Regarding your question about selecting "nfft" value depends on your application/signal. I suggest to study the corresponding chapter in Stoica's spectral estimation book.

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!