Prev Next welch

Spectral Averaging via the Welch-Method
Syntax y = welch(x,Nblk)
Syntax y = welch(x,Nblk,Noffset)
Syntax y = welch(x,Nblk,Noffset,wintype)
Syntax y = welch(x,Nblk,Noffset,wintype,winparam)
Include: include spt\welch.oms
See Also stft

ARGUMENTS:  
   INPUTS:
      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,
                 Noffset<=rowdim(x).
      wintype  = CHAR. One of ["RECTANGULAR"|"TRIANGULAR"|"HAMMING"|"HANNING"|
                 "NUTTALL"|"BLACKMAN"|"BLACKHARRIS3""GAUSSIAN"|"KAISER"].
      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.

Description

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.

Example

Find the spectral density of Additive White Gaussian Noise using Welch's method. Compare this to the un-smoothed spectrum.

include spt\spthead.oms

# 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);