pwelch() - Signal Processing

examplepxx = pwelch(x) returns
the power spectral density (PSD) estimate, pxx,
of the input signal, x, found using Welch's overlapped
segment averaging estimator. When x is a vector,
it is treated as a single channel. When x is
a matrix, the PSD is computed independently for each column and stored
in the corresponding column of pxx. 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.
examplepxx = 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.
examplepxx = pwelch(x,window,noverlap)  uses noverlap samples
of overlap from section to section. noverlap must
be a 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.
examplepxx = 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π).
example[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).

example[___] = 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'.
example[___] = 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'.
example[___] = pwelch(x,window,___,trace) returns
the maximum-hold spectrum estimate if trace is
specified as 'maxhold' and returns the minimum-hold
spectrum estimate if trace is specified as 'minhold'. 

example[___,pxxc] = pwelch(___,'ConfidenceLevel',probability) returns
the probability × 100%
confidence intervals for the PSD estimate in pxxc. 

examplepwelch(___) with no output
arguments plots the Welch PSD estimate in the current figure window.

Syntax

pxx = pwelch(x) examplepxx = pwelch(x,window) examplepxx = pwelch(x,window,noverlap)  examplepxx = pwelch(x,window,noverlap,nfft) example[pxx,w] = pwelch(___)[pxx,f] = pwelch(___,fs) example[pxx,w] = pwelch(x,window,noverlap,w)[pxx,f] = pwelch(x,window,noverlap,f,fs)[___] = pwelch(x,window,___,freqrange) example[___] = pwelch(x,window,___,spectrumtype) example[___] = pwelch(x,window,___,trace) example[___,pxxc] = pwelch(___,'ConfidenceLevel',probability) examplepwelch(___) example

Example

Welch Estimate Using Default InputsOpen This Example
Obtain the Welch PSD estimate of an input signal consisting of a discrete-time sinusoid with an angular frequency of 
 rad/sample with additive 
 white noise.
Create a sine wave with an angular frequency of 
 rad/sample with additive 
 white noise. Reset the random number generator for reproducible results. The signal is 320 samples in length.rng default

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
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 
 rad/sample. Because the signal is real-valued, the periodogram is one-sided and there are 256/2+1 points.pxx = pwelch(x);
plot(10*log10(pxx))

Welch Estimate Using Specified Segment LengthOpen This Example
Obtain the Welch PSD estimate of an input signal consisting of a discrete-time sinusoid with an angular frequency of 
 rad/sample with additive 
 white noise.
Create a sine wave with an angular frequency of 
 rad/sample with additive 
 white noise. Reset the random number generator for reproducible results. The signal is 320 samples in length.rng default

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
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 
 rad/sample. Because the signal is real-valued, the PSD estimate is one-sided and there are 256/2+1 points.segmentLength = 100;
noverlap = 25;
pxx = pwelch(x,segmentLength,noverlap);

plot(10*log10(pxx))

Welch Estimate Specifying Segment OverlapOpen This Example
Obtain the Welch PSD estimate of an input signal consisting of a discrete-time sinusoid with an angular frequency of 
 rad/sample with additive 
 white noise.
Create a sine wave with an angular frequency of 
 rad/sample with additive 
 white noise. Reset the random number generator for reproducible results. The signal is 320 samples in length.rng default

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
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 
 rad/sample. Because the signal is real-valued, the PSD estimate is one-sided and there are 256/2+1 points.segmentLength = 100;
noverlap = 25;
pxx = pwelch(x,segmentLength,noverlap);

plot(10*log10(pxx))

Welch Estimate Using Specified DFT LengthOpen This Example
Obtain the Welch PSD estimate of an input signal consisting of a discrete-time sinusoid with an angular frequency of 
 rad/sample with additive 
 white noise.
Create a sine wave with an angular frequency of 
 rad/sample with additive 
 white noise. Reset the random number generator for reproducible results. The signal is 320 samples in length.rng default

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
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 
 rad/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.segmentLength = 100;
nfft = 640;
pxx = pwelch(x,segmentLength,[],nfft);

plot(10*log10(pxx));
xlabel('Radians/sample')
ylabel('dB')

Welch PSD Estimate of Signal with Frequency in HertzOpen This ExampleCreate a signal consisting of a 100 Hz sinusoid in additive 
 white noise. Reset the random number generator for reproducible results. The sample rate is 1 kHz and the signal is 5 seconds in duration.rng default

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 sample rate to output a vector of frequencies in Hz. Plot the result.[pxx,f] = pwelch(x,500,300,500,fs);

plot(f,10*log10(pxx))

xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')

DC-Centered Power SpectrumOpen This Example
Create a signal consisting of a 100 Hz sinusoid in additive 
 white noise. Reset the random number generator for reproducible results. The sample rate is 1 kHz and the signal is 5 seconds in duration.
rng default

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.[pxx,f] = pwelch(x,500,300,500,fs,'centered','power');

plot(f,10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
grid

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.Maximum-Hold and Minimum-Hold SpectraOpen This Example
Create a signal consisting of three noisy sinusoids and a chirp, sampled at 200 kHz for 0.1 second. The frequencies of the sinusoids are 1 kHz, 10 kHz, and 20 kHz. The sinusoids have different amplitudes and noise levels. The noiseless chirp has a frequency that starts at 20 kHz and increases linearly to 30 kHz during the sampling.
Fs = 200e3;
Fc = [1 10 20]'*1e3;
Ns = 0.1*Fs;

t = (0:Ns-1)/Fs;
x = [1 1/10 10]*sin(2*pi*Fc*t)+[1/200 1/2000 1/20]*randn(3,Ns);
x = x+chirp(t,20e3,t(end),30e3);
Compute the Welch PSD estimate and the maximum-hold and minimum-hold spectra of the signal. Plot the results.[pxx,f] = pwelch(x,[],[],[],Fs);
pmax = pwelch(x,[],[],[],Fs,'maxhold');
pmin = pwelch(x,[],[],[],Fs,'minhold');

plot(f,pow2db(pxx))
hold on
plot(f,pow2db([pmax pmin]),':')
hold off
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
legend('pwelch','maxhold','minhold')

Repeat the procedure, this time computing centered power spectrum estimates.[pxx,f] = pwelch(x,[],[],[],Fs,'centered','power');
pmax = pwelch(x,[],[],[],Fs,'maxhold','centered','power');
pmin = pwelch(x,[],[],[],Fs,'minhold','centered','power');

plot(f,pow2db(pxx))
hold on
plot(f,pow2db([pmax pmin]),':')
hold off
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
legend('pwelch','maxhold','minhold')

Upper and Lower 95%-Confidence BoundsOpen This Example
This 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 
 noise. The amplitude of the two sine waves is 1. The sample rate is 1 kHz. Reset the random number generator for reproducible results.rng default
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-.')

xlim([25 250])
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)');
title('Welch Estimate with 95%-Confidence Bounds');

The lower confidence bound in the immediate vicinity of 100 and 150 Hz is significantly above the upper confidence bound outside the vicinity of 100 and 150 Hz.Welch PSD Estimate of a Multichannel SignalOpen This Example
Generate 1024 samples of a multichannel signal consisting of three sinusoids in additive 
 white Gaussian noise. The sinusoids' frequencies are 
, 
, and 
 rad/sample. Estimate the PSD of the signal using Welch's method and plot it.
N = 1024;
n = 0:N-1;

w = pi./[2;3;4];
x = cos(w*n)' + randn(length(n),3);

pwelch(x)

Related ExamplesBias and Variability in the Periodogram

Output / Return Value

Limitations

Alternatives / See Also

Reference