MATLAB Examples

Streaming Power Spectrum Estimation Using Welch's Method

This example showcases a block that outputs the streaming power spectrum estimate of a time-domain input via Welch's method of averaged modified periodograms.

Contents

Introduction

The spectrum estimator block provides a choice of filter bank-based spectrum estimation and Welch's method of averaged modified periodograms.

In the filter bank method, the input-time doamin signal is divided into different frequency bins using a filter bank and the average power of each subband signal is computed. For more information on filter bank-based spectrum estimation, see High Resolution Spectral Analysis.

In the Welch method, the input time-domain data is partitioned into data segments that are allowed to overlap. A window is applied to each segment, and then an averaged periodogram is computed based on the windowed sequences [1]. The length of the data segments and the choice of the window determine the estimate's resolution bandwidth (RBW), which is the smallest positive frequency that can be resolved.

The block showcased in this example uses the Welch algorithm, which is also used in the Spectrum Analyzer scope block. The Spectrum Analyzer is used to visualize spectra. The block in this example is useful if you need direct access to the estimated spectrum (rather than just visualize it). The output power spectrum may be used as an input to other blocks in your model, or may be logged to the workspace for post-processing.

Block Dialog

This example uses the Spectrum Estimator block.

The block has one input port (the time-domain signal) and up to five output ports:

1. The estimated power spectrum, P.

2. The max-hold spectrum, P_max, which is computed by keeping, at each frequency bin, the maximum value of all the power spectrum estimates. This output is optional.

3. The min-hold spectrum, P_min, which is computed by keeping, at each frequency bin, the minimum value of all the power spectrum estimates. This output is optional.

4. The vector of frequencies, F, at which the power spectrum is estimated. This output is optional.

5. The block's effective resolution bandwidth, RBW. This output is optional.

The block's dialog is shown below. The parameters of this block are equivalent to the ones accessible on the Spectrum Analyzer block. The block dialog allows you to control various aspects of the spectrum estimation algorithm, such as:

  • Select the spectrum type (Power or Power density)
  • Specify how to control the block's resolution (directly specify RBW, or specify a window length).
  • Specify the window function of the estimator.
  • Specify the percentage of overlap between consecutive segments.
  • Select the frequency range as one of two-sided, one-sided, or centered.

Block Architecture

The block is formed of two stages:

1. In the first stage, the input is buffered into segments based on the selected window length and overlap percentage. This stage is implemented using a Buffer block.

2. In the second stage, the power spectrum of the overlapped segments is computed and averaged. This stage is implemented using a dsp.SpectrumEstimator System object.

Specifying Window Length

The dspstreamingwelch model shown below uses a Welch Spectrum Estimator block to estimate the spectrum of a noisy chirp signal sampled at 44100 Hz. The power spectrum estimate is displayed using an Array Plot scope. The peak value of the spectrum, as well as the frequency at which the peak occurs, are detected and displayed on the scope. The estimate's RBW is also displayed. Moreover, a Spectrum Analyzer scope block is included for comparison and validation purposes.

The block's frequency resolution method is set to 'Window Length'. The window length is set to 1024. The FFT length, NFFT, is equal to the window length. The data is windowed using a Chebyshev window with a sidelobe attenuation of 60 dB. The frequency range is one-sided. In this case, the length of the spectrum estimate is NFFT/2+1 = 513 and is computed over the interval [0 Hz,22050 Hz]. The "Sample increment" property of the Array Plot scope is accordingly set to Fs/NFFT = 44100/(1024 * 1000), where the increment is divided by 1000 to scale the frequency units to KHz. You can access the scope's "Sample increment" property by opening its Configuration properties window.

Note that the resolution bandwidth is given by [2] in References:

$$ RBW = enbw(chebwin(N,SL)) * Fs / N $$

where N is the window length, enbw is a function that computes the window's equivalent noise bandwidth, SL is the sidelobe attenuation of the selected Chebyshev window and Fs is the sample rate. In this case, RBW is equal to 65.38 Hz.

When you simulate the model, you can verify that the displayed RBW value is equal to the one shown on the lower bar of the Spectrum Analyzer scope. Moreover, the two blocks give the same peak measurements.

Specifying Non-Zero Overlap

The model in the previous section had zero-overlap. In the dspstreamingwelch_overlap model, we use a Welch Estimation block with an overlap of 50%. Since other model parameters are identical to the previous section, the RBW is unchanged and is equal to 65.38 Hz.

With a the window length of 1024 and an overlap percentage of 50%, 512 input samples are required to form a new data segment. Since the input data is of length 1024, each new data frame yields two new periodograms, and the block's output port runs at a rate twice as fast as the input port.

Note that the Welch estimate block does not have zero latency in this case. The first spectrum estimate output is based on the buffer's initial condition, which is equal to eps. In order to match the spectrum and measurements of the Spectrum Analyzer scope, we therefore insert a delay block at the input of the Spectrum Analyzer.

The results of the Spectrum Analyzer and Welch estimate block may be validated by simulating the model.

Specifying RBW

In the dspstreamingwelch_rbw model, the 'Frequency Resolution Method' parameter is set to RBW. RBW source is Auto. In this mode, similar to the Spectrum Analyzer scope block, the resolution bandwidth is chosen such that there are 1024 RBW intervals over the specified Frequency Span. Since the span in this case is 22050 HZ, the RBW is 21.53 Hz.

The window length used to buffer the data is iteratively computed to yield the desired RBW. The window length in this case is equal to 3073. To verify this value, we can compute the RBW that results from this window length:

$$ RBW = enbw(hann(3073)) * 44100 / 3073 = 21.53 Hz $$

Note that a Hann window is used in this model

In this case, the FFT length, NFFT, is odd and equal to 3073 (the window length). Since the frequency range is one-sided, the spectrum estimate is of length (NFFT + 1)/2 and is computed over the interval [0,44100/2). The "Sample increment" property of the Array Plot scope is set to Fs/NFFT = 44100/(3073 * 1000) KHz.

Again, the results of the Spectrum Analyzer and Welch estimate block may be validated by simulating the model.

References

[1] 'Statistical Digital Signal Processing and Modeling', M.H. Hayes, Wiley, 1996.

[2] http://www.mathworks.com/help/dsp/ref/spectrumanalyzer.html