x = MATRIX, any numerical type. Represents input data by
columns. Coerced to COMPLEX for local processing. Must
have rowdim(x)>=2 or 'novalue' is returned.
Nblk = SCALAR, any numerical type. Block length for spectral
averaging. Coerced to INTEGER for local processing.
Must be Nblk>=2, int(Nblk)<=rowdim(x).
Noffset = SCALAR, any numerical type. Offset for overlapping blocks.
Coerced to INTEGER for local processing. Noffset>=1,
wintype = CHAR. One of ["RECTANGULAR"|"TRIANGULAR"|"HAMMING"|"HANNING"|
winparam = SCALAR, any numerical type. Parameter used by GAUSSIAN or
KAISER window functions. Coerced to DOUBLE for local processing.
RETURN: MATRIX, type DOUBLE, of rowdim=NBLK. Estimated power
spectral density on a column-by-column basis.
Compute Power Spectral Density of an input signal record using Welch's method.
This function estimates the spectral density of an input by segmenting the input record 'x' into blocks of length 'Nblk' samples and averaging their fourier transforms, with optional windowing. The blocks can be overlapped by using the 'Noffset' argument. If the 'Noffset' argument is not used, then the blocks are not overlapped and no windowing is performed. The combination of averaging and overlapping produces a smoother approximation of power density for random signals, though the frequency resolution is more coarse.
Input matrix 'x' represents data arranged as independent columns. The data in fourier transformed into length 'Nblk' blocks, each block starting at integer multiples of 'Noffset', and is processed on a column-by-column basis. Overlapping blocks occur for the case where (Noffset<Nblk), e.g., Noffset=(Nblk/2) implies that later blocks contain the last half of the previous block and contain the first half of the succeeding block.
If the 'wintype'/'winparam' parameters are included, the given window is applied to each block before transformation. Note that applying windows increases the equivalent bandwidth per bin of the resulting spectrum.
Simple row vectors are transposed to column vectors for local processing and returned as row vectors.
The function returns a DOUBLE matrix of row dimension Nblk representing the power spectral density of the input data.
Find the spectral density of Additive White Gaussian Noise using Welch's method. Compare this to the un-smoothed spectrum.
# Form the power spectral desity [PSD] of AWGN
# Use syntax: y = welch(x,Nblk,Noffset, wintype, winparam)
Nrec = 8192; # Samples per input record
fs = 16000d0; # Sampling rate, [Samples/sec]
fc = 1000d0; # Carrier frequency, [Hz]
pc = (33d0)*PI/180d0; # 33 degree Phase offset, in [radians]
x = awgn(0d0, 1d0, Nrec); # Make some noise
# Find the spectral density usng Welch's method
Nblk = 512; # Samples per block
Noffset = 128; # Offset for overlapping (50%)
wintype = "kaiser"; # Choose Kaiser windowing
winparam = 10d0; # Kaiser window factor (2-20)
yWELCH = welch(x,Nblk,Noffset, wintype, winparam); # Find PSD
# Compute the plain FFT of the input record for comparison
yFFT = fft(x);