%% Welch Estimate Using Default Inputs
% 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 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 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.
% Divide the signal into sections of length {\tt nsc}=\lfloor N_{\tt x}/4.5\rfloor. 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=\lceil\log_2{\tt nsc}\rceil.
% 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(:))))

% 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(:))))

% Copyright 2015 The MathWorks, Inc.