Welch’s power spectral density estimate (2024)

Table of Contents

Syntax Description Examples Welch Estimate Using Default Inputs Welch Estimate Using Specified Segment Length Welch Estimate Specifying Segment Overlap Welch Estimate Using Specified DFT Length Welch PSD Estimate of Signal with Frequency in Hertz Maximum-Hold and Minimum-Hold Spectra Upper and Lower 95%-Confidence Bounds DC-Centered Power Spectrum Welch PSD Estimate of a Multichannel Signal Input Arguments noverlap — Number of overlapped samples positive integer | [] nfft — Number of DFT points max(256,2^nextpow2(length(window))) (default) | integer | [] trace — Trace mode 'mean' (default) | 'maxhold' | 'minhold' Output Arguments More About Welch’s Overlapped Segment Averaging Spectral Estimation References Extended Capabilities Tall Arrays Calculate with arrays that have more rows than fit in memory. C/C++ Code Generation Generate C and C++ code using MATLAB® Coder™. Thread-Based EnvironmentRun code in the background using MATLAB® backgroundPool or accelerate code with Parallel Computing Toolbox™ ThreadPool. GPU ArraysAccelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™. Version History R2024a: Code generation support for single-precision variable-size window inputs R2023b: Use single-precision data R2023a: Visualize function outputs using Create Plot Live Editor task See Also Apps Functions Topics MATLAB 命令 Americas Europe Asia Pacific FAQs References

Welch’s power spectral density estimate

collapse all in page

Syntax

pxx = pwelch(x)

pxx = pwelch(x,window)

pxx = pwelch(x,window,noverlap)

pxx = pwelch(x,window,noverlap,nfft)

[pxx,w] = pwelch(___)

[pxx,f] = pwelch(___,fs)

[pxx,w] = pwelch(x,window,noverlap,w)

[pxx,f] = pwelch(x,window,noverlap,f,fs)

[___] = pwelch(x,window,___,freqrange)

[___] = pwelch(x,window,___,trace)

[___,pxxc] = pwelch(___,'ConfidenceLevel',probability)

[___] = pwelch(___,spectrumtype)

pwelch(___)

Description

example

pxx = pwelch(x) returnsthe power spectral density (PSD) estimate, pxx,of the input signal, x, found using Welch's overlappedsegment averaging estimator. When x is a vector,it is treated as a single channel. When x isa matrix, the PSD is computed independently for each column and storedin the corresponding column of pxx. If x isreal-valued, pxx is a one-sided PSD estimate.If x is complex-valued, pxx isa two-sided PSD estimate. By default, x is dividedinto the longest possible segments to obtain as close to but not exceed8 segments with 50% overlap. Each segment is windowed with a Hammingwindow. The modified periodograms are averaged to obtain the PSD estimate.If you cannot divide the length of x exactlyinto an integer number of segments with 50% overlap, x istruncated accordingly.

example

pxx = pwelch(x,window) usesthe input vector or integer, window, to dividethe signal into segments. If window is a vector, pwelch dividesthe signal into segments equal in length to the length of window.The modified periodograms are computed using the signal segments multipliedby the vector, window. If window isan integer, the signal is divided into segments of length window.The modified periodograms are computed using a Hamming window of length window.

example

pxx = pwelch(x,window,noverlap) uses noverlap samplesof overlap from segment to segment. noverlap mustbe a positive integer smaller than window if window isan integer. noverlap must be a positive integerless than the length of window if window isa vector. If you do not specify noverlap, orspecify noverlap as empty, the default numberof overlapped samples is 50% of the window length.

example

pxx = pwelch(x,window,noverlap,nfft) specifiesthe number of discrete Fourier transform (DFT) points to use in thePSD estimate. The default nfft is the greaterof 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 sample rate, 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). fs must be the fifth input to pwelch. To input a sample rate and still use the default values of the preceding optional arguments, specify these arguments as empty, [].

[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 two elements, because otherwise the function interprets it as nfft.

[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 two elements, because otherwise the function interprets it as nfft. The frequencies in f are in cycles per unit time. The sample rate, 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) returnsthe Welch PSD estimate over the frequency range specified by freqrange.Valid options for freqrange are: 'onesided', 'twosided',or 'centered'.

example

[___] = pwelch(x,window,___,trace) returnsthe maximum-hold spectrum estimate if trace isspecified as 'maxhold' and returns the minimum-holdspectrum estimate if trace is specified as 'minhold'.

example

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

example

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

example

pwelch(___) with no outputarguments plots the Welch PSD estimate in the current figure window.

Examples

collapse all

Welch Estimate Using Default Inputs

Open Live Script

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

Create a sine wave with an angular frequency of $π / 4$ rad/sample with additive $N (0, 1)$ white noise. Reset the random number generator for reproducible results. The signal has a length $N_{x} = 320$ samples.

rng defaultn = 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 $2 π / 256$ rad/sample. Because the signal is real-valued, the periodogram is one-sided and there are 256/2+1 points. Plot the Welch PSD estimate.

pxx = pwelch(x);pwelch(x)

Repeat the computation.

Divide the signal into sections of length $n s c = ⌊ N_{x} / 4.5 ⌋$ . This action is equivalent to dividing the signal into the longest possible segments to obtain as close to but not exceed 8 segments with 50% overlap.
Window the sections using a Hamming window.
Specify 50% overlap between contiguous sections
To compute the FFT, use $\max (256, 2^{p})$ points, where $p = ⌈ \log_{2} n s c ⌉$ .

Verify that the two approaches give identical results.

Nx = length(x);nsc = floor(Nx/4.5);nov = floor(nsc/2);nff = max(256,2^nextpow2(nsc));t = pwelch(x,hamming(nsc),nov,nff);maxerr = max(abs(abs(t(:))-abs(pxx(:))))

maxerr = 0

Divide the signal into 8 sections of equal length, with 50% overlap between sections. Specify the same FFT length as in the preceding step. Compute the Welch PSD estimate and verify that it gives the same result as the previous two procedures.

ns = 8;ov = 0.5;lsc = floor(Nx/(ns-(ns-1)*ov));t = pwelch(x,lsc,floor(ov*lsc),nff);maxerr = max(abs(abs(t(:))-abs(pxx(:))))

maxerr = 0

Welch Estimate Using Specified Segment Length

Open Live Script

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

Create a sine wave with an angular frequency of $π / 3$ rad/sample with additive $N (0, 1)$ white noise. Reset the random number generator for reproducible results. The signal has 512 samples.

rng defaultn = 0:511;x = cos(pi/3*n)+randn(size(n));

Obtain the Welch PSD estimate dividing the signal into segments 132 samples in length. The signal segments are multiplied by a Hamming window 132 samples in length. The number of overlapped samples is not specified, so it is set to 132/2 = 66. The DFT length is 256 points, yielding a frequency resolution of $2 π / 256$ rad/sample. Because the signal is real-valued, the PSD estimate is one-sided and there are 256/2+1 = 129 points. Plot the PSD as a function of normalized frequency.

segmentLength = 132;[pxx,w] = pwelch(x,segmentLength);plot(w/pi,10*log10(pxx))xlabel('\omega / \pi')

Welch Estimate Specifying Segment Overlap

Open Live Script

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

Create a sine wave with an angular frequency of $π / 4$ rad/sample with additive $N (0, 1)$ white noise. Reset the random number generator for reproducible results. The signal is 320 samples in length.

rng defaultn = 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 $2 π / 256$ 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 Length

Open Live Script

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

rng defaultn = 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 $π / 4$ 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('rad/sample')ylabel('dB / (rad/sample)')

Welch PSD Estimate of Signal with Frequency in Hertz

Open Live Script

Create a signal consisting of a 100 Hz sinusoid in additive N(0,1) 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 defaultfs = 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('PSD (dB/Hz)')

Maximum-Hold and Minimum-Hold Spectra

Open Live Script

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 onplot(f,pow2db([pmax pmin]),':')hold offxlabel('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 onplot(f,pow2db([pmax pmin]),':')hold offxlabel('Frequency (Hz)')ylabel('Power (dB)')legend('pwelch','maxhold','minhold')

Upper and Lower 95%-Confidence Bounds

Open Live Script

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 N(0,1) 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 defaultfs = 1000;t = 0:1/fs:1-1/fs;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 onplot(f,10*log10(pxxc),'-.')hold offxlim([25 250])xlabel('Frequency (Hz)')ylabel('PSD (dB/Hz)')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.

DC-Centered Power Spectrum

Open Live Script

Create a signal consisting of a 100 Hz sinusoid in additive $N (0, 1 / 4)$ 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 defaultfs = 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.

Welch PSD Estimate of a Multichannel Signal

Open Live Script

Generate 1024 samples of a multichannel signal consisting of three sinusoids in additive $N (0, 1)$ white Gaussian noise. The sinusoids' frequencies are $π / 2$ , $π / 3$ , and $π / 4$ 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)

Input Arguments

collapse all

`noverlap` — Number of overlapped samples
positive integer | `[]`

Number of overlapped samples, specified as a positive integersmaller than the length of window. If you omit noverlap orspecify noverlap as empty, a value is used toobtain 50% overlap between segments.

`nfft` — Number of DFT points
`max(256,2^nextpow2(length(window)))` (default) | integer | `[]`

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

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

Data Types: single | double

`trace` — Trace mode
`'mean'` (default) | `'maxhold'` | `'minhold'`

Trace mode, specified as one of 'mean', 'maxhold',or 'minhold'. The default is 'mean'.

'mean' — returns the Welchspectrum estimate of each input channel. pwelch computesthe Welch spectrum estimate at each frequency bin by averaging thepower spectrum estimates of all the segments.
'maxhold' — returns themaximum-hold spectrum of each input channel. pwelch computesthe maximum-hold spectrum at each frequency bin by keeping the maximumvalue among the power spectrum estimates of all the segments.
'minhold' — returns theminimum-hold spectrum of each input channel. pwelch computesthe minimum-hold spectrum at each frequency bin by keeping the minimumvalue among the power spectrum estimates of all the segments.

Output Arguments

collapse all

More About

collapse all

Welch’s Overlapped Segment Averaging Spectral Estimation

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.

References

[1] Hayes, Monson H. Statistical Digital Signal Processing and Modeling. New York: John Wiley & Sons, 1996.

[2] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.

Extended Capabilities

Tall Arrays
Calculate with arrays that have more rows than fit in memory.

Usage notes and limitations:

The input x must not be a tall row vector
The window argument must always be specified.

For more information, see Tall Arrays.

C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.

Thread-Based Environment
Run code in the background using MATLAB® `backgroundPool` or accelerate code with Parallel Computing Toolbox™ `ThreadPool`.

Usage notes and limitations:

The syntax with no output arguments is not supported.

For more information, see Run MATLAB Functions in Thread-Based Environment.

GPU Arrays
Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

This function fully supports GPU arrays. For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced before R2006a

expand all

R2024a: Code generation support for single-precision variable-size window inputs

The pwelch function supports single-precision variable-size window inputs for code generation.

R2023b: Use single-precision data

The pwelch function supports single-precision inputs.

R2023a: Visualize function outputs using Create Plot Live Editor task

You can now use the Create Plot Live Editor task to visualize the output of pwelch interactively. You can select different chart types and set optional parameters. The task also automatically generates code that becomes part of your live script.

MATLAB 命令

您点击的链接对应于以下 MATLAB 命令：

请在 MATLAB 命令行窗口中直接输入以执行命令。Web 浏览器不支持 MATLAB 命令。

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

Americas

América Latina (Español)
Canada (English)
United States (English)

Europe

Belgium (English)
Denmark (English)
Deutschland (Deutsch)
España (Español)
Finland (English)
France (Français)
Ireland (English)
Italia (Italiano)
Luxembourg (English)

Netherlands (English)
Norway (English)
Österreich (Deutsch)
Portugal (English)
Sweden (English)
Switzerland
United Kingdom (English)

Asia Pacific

Australia (English)
India (English)
New Zealand (English)
中国
- 简体中文
- English
日本 (日本語)
한국 (한국어)

Contact your local office

Welch’s power spectral density estimate (2024)

FAQs

What is Welch's power spectral density estimate? ›

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.

Tell Me More ›

What is the Welch method in EEG? ›

The Welch method is another commonly used to estimate PSD from EEG. The basic idea behind the Welch method is to use a moving window technique where the FFT is computed in each window and the PSD is then computed as an average of FFTs over all windows. This is achieved by using the pwelch() command in MATLAB.

Read On ›

What is the Welch method in DSP? ›

The Welch method reduces the variance of the periodogram method by averaging. This method first divides a time series into overlapping subsequences by applying a window to each subsequence and then averaging the periodogram of each subsequence.

Get More Info ›

What does power spectral density tell us? ›

A Power Spectral Density (PSD) is the measure of signal's power content versus frequency. A PSD is typically used to characterize broadband random signals. The amplitude of the PSD is normalized by the spectral resolution employed to digitize the signal.

Get More Info Here ›

When to use the Welch method? ›

It is used in physics, engineering, and applied mathematics for estimating the power of a signal at different frequencies. The method is based on the concept of using periodogram spectrum estimates, which are the result of converting a signal from the time domain to the frequency domain.

What is the difference between power spectral density and FFT? ›

The PSD and FFT are tools for measuring and analyzing a signal's frequency content. The FFT transfers time data to the frequency domain, which allows engineers to view changes in frequency values. The PSD takes another step and calculates the power, or strength, of the frequency content.

Learn More ›

What is the interpretation of spectral density? ›

Intuitively speaking, the spectral density characterizes the frequency content of the signal. One purpose of estimating the spectral density is to detect any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities.

Tell Me More ›

What is the problem of spectral estimation? ›

Spectral estimation is the problem of estimating the power spectrum of a stochastic process given partial data, usually only a finite number of samples of the autocorrelation function of limited accuracy.

What is Welch approximation? ›

In statistics and uncertainty analysis, the Welch–Satterthwaite equation is used to calculate an approximation to the effective degrees of freedom of a linear combination of independent sample variances, also known as the pooled degrees of freedom, corresponding to the pooled variance.

Learn More ›

What is the formula for power spectral density? ›

4 Power Spectral Density and Dynamic Range. A signal consisting of many similar subcarriers will have a constant power spectral density (PSD) over its bandwidth and the total signal power can then be found as P = PSD · BW.

Discover More Details ›

What is the difference between pwelch and periodogram? ›

In periodogram, you sample the [0 Fs] with 1000 points, so the bin width, i.e., the frequency resolution for each sampling point, is 1 Hz. In pwelch, since you are using only 500 points, each bin becomes 2 Hz.

Discover More Details ›

What is Welch's method of averaged modified periodogram? ›

Welch's method [1] computes an estimate of the power spectral density by dividing the data into overlapping segments, computing a modified periodogram for each segment and averaging the periodograms. Sampling frequency of the x time series.