How to use pwelch for PSD?

12 visualizaciones (últimos 30 días)
Louise Wilson
Louise Wilson el 26 de Jul. de 2019
Respondida: emehmetcik el 8 de Nov. de 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!

Respuestas (1)

emehmetcik
emehmetcik el 8 de Nov. de 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.

Categorías

Más información sobre Parametric Spectral Estimation en Help Center y File Exchange.

Productos


Versión

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by