Welch’s power spectral density estimate (2024)

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

Create a sine wave with an angular frequency of $$\pi /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\pi /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.

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

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

Create a sine wave with an angular frequency of $$\pi /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\pi /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')

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

Create a sine wave with an angular frequency of $$\pi /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\pi /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))

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

Create a sine wave with an angular frequency of $$\pi /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. Use the default overlap of 50%. Specify the DFT length to be 640 points so that the frequency of $$\pi /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)')

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)')

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')

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.

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.

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