pwelch (Signal Processing Toolbox)

Estimate the power spectral density (PSD) of a signal using Welch's method.

Syntax

Pxx = pwelch(x)
Pxx = pwelch(x,nfft)
[Pxx,w] = pwelch(x,nfft)
[Pxx,f] = pwelch(x,nfft,Fs)
[Pxx,f] = pwelch(x,nfft,Fs,window)
[Pxx,f] = pwelch(x,nfft,Fs,window,noverlap)
[Pxx,Pxxc,f] = pwelch(x,nfft,Fs,window,noverlap,p)
[Pxx,Pxxc,f] = pwelch(...,'range')
pwelch(...)
pwelch(...,'magunits')

Description

Pxx = pwelch(x)

estimates the power spectrum of the sequence x using the Welch method of spectral estimation. If x is real, pwelch estimates the spectrum at positive frequencies only; in this case, output Pxx is a column vector of length nfft/2+1 for nfft even and (nfft+1)/2 for nfft odd. If x is complex, pwelch estimates the spectrum at both positive and negative frequencies and Pxx has length nfft.

Pxx = pwelch(x,nfft)

uses the specified FFT length nfft in estimating the power spectrum for x. This value determines the number of different frequencies at which the power spectrum is estimated. Specify nfft as a power of 2 for fastest execution. Specify an empty matrix for nfft, [], to use the default value of min(256,length(x)).

[Pxx,w] = pwelch(x,nfft)

returns a vector w of normalized angular frequencies (in rads/sample) at which the function evaluates the PSD. The range of w is [0,

] for real x and [0,2

] for complex x. Since the frequency vector w is the same size as Pxx, plot(w,Pxx) plots the power spectrum versus the normalized angular frequency.

[Pxx,f] = pwelch(x,nfft,Fs)

returns a vector f of linear frequencies (in Hz) at which the function evaluates the PSD. Fs is a scalar that specifies the sampling frequency of x, and Pxx is scaled by 1/Fs. The range of f is [0,Fs/2] for real x and [0,Fs] for complex x. Since the frequency vector f is the same size as Pxx, plot(f,Pxx) plots the power spectrum versus the linear frequency. Specify an empty matrix for Fs, [], to use the default value of 1 Hz.

[Pxx,f] = pwelch(x,nfft,Fs,window)

specifies a windowing function and the number of samples per windowed section of the x vector (e.g., kaiser(64)). The length of the window must be less than or equal to nfft; pwelch zero pads the sections if the length of the window is less than nfft. Specify a scalar for window to use a Hanning window of that length; specify an empty matrix, [], to use the default value of hanning(nfft).

[Pxx,f] = pwelch(x,nfft,Fs,window,noverlap)

overlaps the windowed sections of x by noverlap samples. Specify an empty matrix for noverlap, [], to use the default value of 0.

[Pxx,Pxxc,f] = pwelch(x,nfft,Fs,window,noverlap,p)

where p is a positive scalar between 0 and 1 returns a vector Pxxc that contains an estimate of the p*100 percent confidence interval for Pxx. Pxxc is a two-column matrix that is the same length as Pxx. The interval [Pxxc(:,1),Pxxc(:,2)] covers the true PSD with probability p. plot(f,[Pxx Pxxc]) plots the power spectrum inside the p*100 percent confidence interval. Specify an empty matrix for p, [], to use the default value of 0.95.

[Pxx,Pxxc,f] = pwelch(...,'range')

specifies the range of frequency values to include in f. range can be:

half, to compute the PSD over the range [0,] for real or complex x. If Fs is specified, the range is [0,Fs/2] for real or complex x. If Fs is specified by an empty matrix, [], the default value of 1 is used, and the range is [0,1/2] for real or complex x.
whole, to compute the PSD over the range [0,2] for real or complex x. If Fs is specified, the range is [0,Fs] for real or complex x. If Fs is specified by an empty matrix, [], the default value of 1 is used, and the range is [0,1] for real or complex x.

pwelch(x,...)

with no output arguments plots the PSD against frequency in the current figure window. If p is specified, the plot includes the confidence interval.

pwelch(...,'magunits')

specifies the units in which to express to the magnitude axis of the PSD plot. The options are:

'db', to express the PSD data in decibels (dB). This is the default.
'squared', to express the PSD in linear units (magnitude squared).

Example

Generate a colored noise signal and plot its PSD with a confidence interval of 95%. Specify a length 1024 FFT, a 512-point Kaiser window with no overlap, and a sampling frequency of 10 kHz:

h = fir1(30,0.2,boxcar(31));     % design a lowpass filter
r = randn(16384,1);              % white noise
x = filter(h,1,r);               % color the noise
pwelch(x,1024,10000,kaiser(512,5),0,0.95)

Algorithm

pwelch calculates the power spectral density using Welch's method (see references [1] and [2]):

1.: It applies the window specified by the window vector to each successive section of input x.

2.: It transforms each section with an nfft-point FFT.

3.: It forms the periodogram of each section by scaling the magnitude squared of each transformed section.

4.: It averages the periodograms of the overlapping sections to form Pxx(f).

5.: It scales Pxx(f) by 1/Fs to form Pxx(f)/Fs, the power spectrum of x.

The number of sections that pwelch averages is

k = fix((length(x)-noverlap)/(length(window)-noverlap))

Diagnostics

An appropriate diagnostic message is displayed when incorrect arguments to pwelch are used:

Requires window's length to be no greater than FFT length.
Requires NOVERLAP to be strictly less than the window length.
Requires positive integer values for NFFT and NOVERLAP.
Requires confidence parameter to be a scalar between 0 and 1.
Requires vector input.

`cohere`	Estimate magnitude squared coherence function between two signals.
`csd`	Estimate the cross spectral density (CSD) of two signals.
`pburg`	Power spectrum estimate using the Burg method.
`pcov`	Power spectrum estimate using the covariance method.
`pmcov`	Power spectrum estimate using the modified covariance method.
`pmtm`	Power spectrum estimate using the multitaper method (MTM).
`pmusic`	Power spectrum estimate using MUSIC eigenvector method.
`pyulear`	Power spectrum estimate using Yule-Walker AR method.
`specgram`	Time-dependent frequency analysis (spectrogram).
`tfe`	Transfer function estimate from input and output.