pxx = pwelch(x) returns
the power spectral density (PSD) estimate, pxx,
of the input signal, x using Welch's overlapped
segment averaging estimator. If x is real-valued, pxx is
a one-sided PSD estimate. If x is complex-valued, pxx is
a two-sided PSD estimate. By default, x is divided
into the longest possible sections to obtain as close to but not exceed
8 segments with 50% overlap. Each section is windowed with a Hamming
window. The modified periodograms are averaged to obtain the PSD estimate.
If you cannot divide the length of x exactly
into an integer number of sections with 50% overlap, x is
truncated accordingly.

pxx = pwelch(x,window) uses
the input vector or integer, window, to divide
the signal into sections. If window is a vector, pwelch divides
the signal into sections equal in length to the length of window.
The modified periodograms are computed using the signal sections multiplied
by the vector, window. If window is
an integer, the signal is divided into sections of length window.
The modified periodograms are computed using a Hamming window of length window.

pxx = pwelch(x,window,noverlap) uses noverlap samples
of overlap from section to section. noverlap must
be an positive integer smaller than window if window is
an integer. noverlap must be a positive integer
less than the length of window if window is
a vector. If you do not specify noverlap, or
specify noverlap as empty, the default number
of overlapped samples is 50% of the window length.

pxx = pwelch(x,window,noverlap,nfft) specifies
the number of discrete Fourier transform (DFT) points to use in the
PSD estimate. The default nfft is the greater
of 256 or the next power of 2 greater than the length of the segments.

[pxx,w]
= pwelch(___)returns the normalized frequency
vector, w. If pxx is a one-sided
PSD estimate, w spans the interval [0,π]
if nfft is even and [0,π) if nfft is
odd. If pxx is a two-sided PSD estimate, w spans
the interval [0,2π).

[pxx,f]
= pwelch(___,fs) returns
a frequency vector, f, in cycles per unit time.
The sampling frequency, fs, is the number of
samples per unit time. If the unit of time is seconds, then f is
in cycles/sec (Hz). For real–valued signals, f spans
the interval [0,fs/2] when nfft is
even and [0,fs/2) when nfft is
odd. For complex-valued signals, f spans the
interval [0,fs).

[pxx,w]
= pwelch(x,window,noverlap,w) returns
the two-sided Welch PSD estimates at the normalized frequencies specified
in the vector, w. The vector, w,
must contain at least 2 elements.

[pxx,f]
= pwelch(x,window,noverlap,f,fs) returns
the two-sided Welch PSD estimates at the frequencies specified in
the vector, f. The vector, f,
must contain at least 2 elements. The frequencies in f are
in cycles per unit time. The sampling frequency, fs,
is the number of samples per unit time. If the unit of time is seconds,
then f is in cycles/sec (Hz).

[___] = pwelch(x,window,___,freqrange) returns
the Welch PSD estimate over the frequency range specified by freqrange.
Valid options for freqrange are: 'onesided', 'twosided',
or 'centered'.

[___] = pwelch(x,window,___,spectrumtype) returns
the PSD estimate if spectrumtype is specified
as 'psd' and returns the power spectrum if spectrumtype is
specified as 'power'.

Obtain the Welch PSD estimate of an input signal
consisting of a discrete-time sinusoid with an angular frequency of
π/4 radians/sample with additive N(0,1) white noise.

Create a sine wave with an angular frequency of π/4
radians/sample with additive N(0,1) white noise. The signal is 320
samples in length. Obtain the Welch PSD estimate using the default
Hamming window and DFT length. The default segment length is 71 samples
and the DFT length is the 256 points yielding a frequency resolution
of 2π/256 radians/sample. Because the signal is real-valued,
the periodogram is one-sided and there are 256/2+1 points.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
pxx = pwelch(x);
plot(10*log10(pxx))

Obtain the Welch PSD estimate of an input signal
consisting of a discrete-time sinusoid with an angular frequency of
π/4 radians/sample with additive N(0,1) white noise.

Create a sine wave with an angular frequency of π/4
radians/sample with additive N(0,1) white noise. The signal is 320
samples in length. Obtain the Welch PSD estimate dividing the signal
into segments 100 samples in length. The signal segments are multiplied
by a Hamming window 100 samples in length. The number of overlapped
samples is 50. The DFT length is 256 points yielding a frequency resolution
of 2π/256 radians/sample. Because the signal is real-valued,
the PSD estimate is one-sided and there are 256/2+1 points.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
segmentLength = 100;
pxx = pwelch(x,segmentLength);
plot(10*log10(pxx))

Obtain the Welch PSD estimate of an input signal
consisting of a discrete-time sinusoid with an angular frequency of
π/4 radians/sample with additive N(0,1) white noise.

Create a sine wave with an angular frequency of π/4
radians/sample with additive N(0,1) white noise. The signal is 320
samples in length. Obtain the Welch PSD estimate dividing the signal
into segments 100 samples in length. The signal segments are multiplied
by a Hamming window 100 samples in length. The number of overlapped
samples is 25. The DFT length is 256 points yielding a frequency resolution
of 2π/256 radians/sample. Because the signal is real-valued,
the PSD estimate is one-sided and there are 256/2+1 points.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
segmentLength = 100;
noverlap = 25;
pxx = pwelch(x,segmentLength,noverlap);
plot(10*log10(pxx))

Obtain the Welch PSD estimate of an input signal
consisting of a discrete-time sinusoid with an angular frequency of
π/4 radians/sample with additive N(0,1) white noise.

Create a sine wave with an angular frequency of π/4
radians/sample with additive N(0,1) white noise. The signal is 320
samples in length. Obtain the Welch PSD estimate dividing the signal
into segments 100 samples in length. Use the default overlap of 50%.
Specify the DFT length to be 640 points so that the frequency of π/4
radians/sample corresponds to a DFT bin (bin 81).Because the signal
is real-valued, the PSD estimate is one-sided and there are 640/2+1
points.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
segmentLength = 100;
nfft = 640;
pxx = pwelch(x,segmentLength,[],nfft);
plot(10*log10(pxx));
xlabel('Radians/sample'); ylabel('dB');

Create a signal consisting of a 100-Hz sinusoid in additive
N(0,1) white noise. The sampling rate is 1 kHz and the signal is 5
seconds in duration.

fs = 1000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t)+randn(size(t));

Obtain Welch's overlapped segment averaging PSD
estimate of the preceding signal. Use a segment length of 500 samples
with 300 overlapped samples. Use 500 DFT points so that 100 Hz falls
directly on a DFT bin. Input the sampling frequency to output a vector
of frequencies in Hz. Plot the result.

Create a signal consisting of a 100-Hz sinusoid in additive
N(0,1/4) white noise. The sampling rate is 1 kHz and the signal is
5 seconds in duration.

fs = 1000;
t = 0:1/fs:5-1/fs;
noisevar = 1/4;
x = cos(2*pi*100*t)+sqrt(noisevar)*randn(size(t));

Obtain the DC-centered power spectrum using Welch's
method. Use a segment length of 500 samples with 300 overlapped samples
and a DFT length of 500 points. Plot the result.

You see that the power at –100 and 100 Hz is close to
the expected power of 1/4 for a real-valued sine wave with an amplitude
of 1. The deviation from 1/4 is due to the effect of the additive
noise.

The following example illustrates the use of
confidence bounds with Welch's overlapped segment averaging
(WOSA) PSD estimate. While not a necessary condition for statistical
significance, frequencies in Welch's estimate where the lower
confidence bound exceeds the upper confidence bound for surrounding
PSD estimates clearly indicate significant oscillations in the time
series.

Create a signal consisting of the superposition of 100-Hz
and 150-Hz sine waves in additive white N(0,1) noise. The amplitude
of the two sine waves is 1. The sampling frequency is 1 kHz.

t = 0:0.001:1-0.001;
fs = 1000;
x = cos(2*pi*100*t)+sin(2*pi*150*t)+randn(size(t));

Obtain the WOSA estimate with 95%-confidence bounds. Set
the segment length equal to 200 and overlap the segments by 50% (100
samples). Plot the WOSA PSD estimate along with the confidence interval
and zoom in on the frequency region of interest near 100 and 150 Hz.

L = 200;
noverlap = 100;
[pxx,f,pxxc] = pwelch(x,hamming(L),noverlap,200,fs,...'ConfidenceLevel',0.95);
plot(f,10*log10(pxx)); hold on;
plot(f,10*log10(pxxc),'r--','linewidth',2);
axis([25 250 min(min(10*log10(pxxc))) max(max(10*log10(pxxc)))]);
xlabel('Hz'); ylabel('dB');
title('Welch Estimate with 95%-Confidence Bounds');

At 100 and 150 Hz, the lower confidence bound exceeds the upper
confidence bounds for surrounding PSD estimates.

Window, specified as a row or column vector or an integer. If window is
a vector, pwelch divides x into
overlapping sections of length equal to the length of window,
and then multiplies each signal section with the vector specified
in window. If window is
an integer, pwelch is divided into sections of
length equal to the integer value, and a Hamming window of equal length
is used. If the length of
x cannot be divided exactly into an integer number
of sections with noverlap number of overlapping
samples, x is truncated accordingly. If you specify window as
empty, the default Hamming window is used to obtain eight sections
of x with noverlap overlapping
samples.

Number of overlapped samples, specified as a positive integer
smaller than the length of window. If you omit noverlap or
specify noverlap as empty, a value is used to
obtain 50% overlap between segments.

Number of DFT points, specified as a positive integer. For a
real-valued input signal, x, the PSD estimate, pxx has
length (nfft/2+1) if nfft is
even, and (nfft+1)/2 if nfft is
odd. For a complex-valued input signal,x, the
PSD estimate always has length nfft. If nfft is
specified as empty, the default nfft is used.

If nfft is greater than the segment length,
the data is zero-padded. If nfft is less than
the segment length, the segment is wrapped using datawrap to
make the length equal to nfft.

Sampling frequency specified as a positive scalar. The sampling
frequency is the number of samples per unit time. If the unit of time
is seconds, the sampling frequency has units of hertz.

Normalized frequencies for Goertzel algorithm, specified as
a row or column vector with at least 2 elements. Normalized frequencies
are in radians/sample.

Cyclical frequencies for Goertzel algorithm, specified as a
row or column vector with at least 2 elements. The frequencies are
in cycles per unit time. The unit time is specified by the sampling
frequency, fs. If fs has
units of samples/second, then f has units of
Hz.

Frequency range for the PSD estimate, specified as a one of 'onesided', 'twosided',
or 'centered'. The default is 'onesided' for
real-valued signals and 'twosided' for complex-valued
signals. The frequency ranges corresponding to each option are

'onesided' — returns the
one-sided PSD estimate of a real-valued input signal, x.
If nfft is even, pxx will
have length nfft/2+1 and is computed over the
interval [0,π] radians/sample. If nfft is
odd, the length of pxx is (nfft+1)/2
and the interval is [0,π) radians/sample. When fs is
optionally specified, the corresponding intervals are [0,fs/2]
cycles/unit time and [0,fs/2) cycles/unit time
for even and odd length nfft respectively.

'twosided' — returns the
two-sided PSD estimate for either the real-valued or complex-valued
input, x. In this case, pxx has
length nfft and is computed over the interval
[0,2π) radians/sample. When fs is optionally
specified, the interval is [0,fs) cycles/unit
time.

'centered' — returns the
centered two-sided PSD estimate for either the real-valued or complex-valued
input, x. In this case, pxx has
length nfft and is computed over the interval
(-π,π] radians/sample for even length nfft and
(-π,π) radians/sample for odd length nfft.
When fs is optionally specified, the corresponding
intervals are (-fs/2, fs/2]
cycles/unit time and (-fs/2, fs/2)
cycles/unit time for even and odd length nfft respectively.

Power spectrum scaling, specified as one of 'psd' or 'power'.
Omitting the spectrumtype, or specifying 'psd',
returns the power spectral density. Specifying 'power' scales
each estimate of the PSD by the equivalent noise bandwidth of the
window. Use the 'power' option to obtain an estimate
of the power at each frequency.

Coverage probability for the true PSD, specified as a scalar
in the range (0,1). The output, pxxc, contains
the lower and upper bounds of the probabilityx100%
interval estimate for the true PSD.

PSD estimate, specified as a real-valued, nonnegative column
vector. The units of the PSD estimate are in squared magnitude units
of the time series data per unit frequency. For example, if the input
data is in volts, the PSD estimate is in units of squared volts per
unit frequency. For a time series in volts, if you assume a resistance
of 1 ohm and specify the sampling frequency in Hz, the PSD estimate
is in watts/Hz.

Normalized frequencies, specified as a real-valued column vector.
If pxx is a one-sided PSD estimate, w spans
the interval [0,π] if nfft is even and
[0,π) if nfft is odd. If pxx is
a two-sided PSD estimate, w spans the interval
[0,2π). For a DC-centered PSD estimate, f spans
the interval (-π,π] radians/sample for even length nfft and
(-π,π) radians/sample for odd length nfft.

Cyclical frequencies, specified as a real-valued column vector.
For a one-sided PSD estimate, f spans the interval
[0,fs/2] when nfft is even
and [0,fs/2) when nfft is
odd. For a two-sided PSD estimate, f spans the
interval [0,fs). For a DC-centered PSD estimate, f spans
the interval (-fs/2, fs/2]
cycles/unit time for even length nfft and (-fs/2, fs/2)
cycles/unit time for odd length nfft .

Confidence bounds, specified as an N-by-2 matrix with real-valued
elements. The row dimension of the matrix is equal to the length of
the PSD estimate, pxx. The first column contains
the lower confidence bound and the second column contains the upper
confidence bound for the corresponding PSD estimates in the rows of pxx.
The coverage probability of the confidence intervals is determined
by the value of the probability input.

The periodogram is not a consistent estimator of the true power
spectral density of a wide-sense stationary process. Welch's
technique to reduce the variance of the periodogram breaks the time
series into segments, usually overlapping. Welch's method computes
a modified periodogram for each segment and then averages these estimates
to produce the estimate of the power spectral density. Because the
process is wide-sense stationary and Welch's method uses PSD
estimates of different segments of the time series, the modified periodograms
represent approximately uncorrelated estimates of the true PSD and
averaging reduces the variability.

The segments are typically multiplied by a window function,
such as a Hamming window, so that Welch's method amounts to
averaging modified periodograms. Because the segments usually overlap,
data values at the beginning and end of the segment tapered by the
window in one segment, occur away from the ends of adjacent segments.
This guards against the loss of information caused by windowing.