dsp.SpectrumEstimator

Estimate power spectrum or power density spectrum

Description

The dsp.SpectrumEstimator System object™ computes the power spectrum or the power density spectrum of a signal using the Welch algorithm or the filter bank approach.

When you choose the Welch method, the object computes the averaged modified periodograms to compute the spectral estimate. When you choose the filter bank approach, an analysis filter bank splits the broadband input signal into multiple narrow subbands. The object computes the power in each narrow frequency band, and the computed value is the spectral estimate over the respective frequency band. For signals with relatively small FFT lengths, the filter bank approach produces a spectral estimate with a higher resolution, a more accurate noise floor, and peaks more precise than the Welch method, with low or no spectral leakage. These advantages come at the expense of increased computation and slower tracking.

The spectrum can be expressed in watts or in decibels. This object can also estimate the max-hold and min-hold spectra of the signal.

To estimate the power density spectrum:

1. Create the dsp.SpectrumEstimator object and set its properties.

2. Call the object with arguments, as if it were a function.

To learn more about how System objects work, see What Are System Objects?

Creation

Syntax

SE = dsp.SpectrumEstimator

SE = dsp.SpectrumEstimator(Name,Value)

Description

SE = dsp.SpectrumEstimator returns a System object, SE, that computes the frequency power spectrum or the power density spectrum of real or complex signals. This System object uses the Welch’s averaged modified periodogram method or the filter bank-based spectral estimation method.

SE = dsp.SpectrumEstimator(Name,Value) returns a dsp.SpectrumEstimator System object with each specified property name set to the specified value. Unspecified properties have default values.

Properties

Unless otherwise indicated, properties are nontunable, which means you cannot change their values after calling the object. Objects lock when you call them, and the release function unlocks them.

If a property is tunable, you can change its value at any time.

For more information on changing property values, see System Design in MATLAB Using System Objects.

SpectrumType — Spectrum type

'Power' (default) | 'Power density'

Spectrum type, specified as either 'Power' or 'Power density'. When the spectrum type is 'Power', the power density spectrum is scaled by the equivalent noise bandwidth of the window (in Hz).

Tunable: Yes

FFTLengthSource — Source of FFT length value

'Auto' (default) | 'Property'

Source of the FFT length value, specified as either 'Auto' or 'Property'. If you set this property to 'Auto', the spectrum estimator sets the FFT length to the input frame size. If you set this property to 'Property', then you specify the number of FFT points using the FFTLength property.

FFTLength — FFT length

128 (default) | positive integer

Specify the length of the FFT that the spectrum estimator uses to compute spectral estimates as a positive integer.

Dependencies

This property applies when you set the FFTLengthSource property to 'Property'.

Data Types: double

Method — Welch or filter bank

'Welch' (default) | 'Filter bank'

Specify the spectral estimation method:

• 'Welch' — The object uses Welch's averaged modified periodograms method.

• 'Filter bank' — An analysis filter bank splits the broadband input signal into multiple narrow subbands. The object computes the power in each narrow frequency band, and the computed value is the spectral estimate over the respective frequency band.

Window — Window function

'Hann' (default) | 'Rectangular' | 'Chebyshev' | 'Flat Top' | 'Hamming' | 'Kaiser'

Specify a window function for the spectral estimator as one of 'Rectangular', 'Chebyshev', 'Flat Top', 'Hamming', 'Hann', or 'Kaiser'.

Dependencies

This property applies when you set Method to 'Welch'.

NumTapsPerBand — Number of filter taps per frequency band

12 (default) | positive integer

Specify the number of filter coefficients, or taps, for each frequency band. This value corresponds to the number of filter coefficients per polyphase branch. The total number of filter coefficients is given by NumTapsPerBand × FFTLength.

Dependencies

This property applies when you set Method to 'Filter bank'.

Data Types: single | double | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64

FrequencyRange — Frequency range of the spectrum estimate

'twosided' (default) | 'onesided' | 'centered'

Specify the frequency range of the spectrum estimator as one of 'twosided', 'onesided', or 'centered'.

If you set the FrequencyRange to 'onesided', then the spectrum estimator computes the one-sided spectrum of a real input signal. When the FFT length, NFFT, is even, the spectrum estimate has length (NFFT/2) + 1 and is computed over the frequency range [0, SampleRate/2], where SampleRate is the sample rate of the input signal. When NFFT is odd, the spectrum estimate has length (NFFT + 1)/2 and is computed over the frequency range [0, SampleRate/2).

If you set the FrequencyRange to 'twosided', then the spectrum estimator computes the two-sided spectrum of a complex or real input signal. The length of the spectrum estimate is equal to the FFT length. The spectrum estimate is computed over the frequency range [0, SampleRate), where SampleRate is the sample rate of the input signal.

If you set the FrequencyRange to 'centered', then the spectrum estimator computes the centered two-sided spectrum of a complex or real input signal. The length of the spectrum estimate is equal to the FFT length. The spectrum estimate is computed over the frequency range (-SampleRate/2, SampleRate/2] when the FFT length is even and (-SampleRate/2, SampleRate/2) when the FFT length is odd.

PowerUnits — Power units

'Watts' (default) | 'dBW' | 'dBm'

Specify the units used to measure power as one of 'Watts', 'dBW', or 'dBm'.

AveragingMethod — Averaging method

'Running' (default) | 'Exponential'

Specify the averaging method as 'Running' or 'Exponential'. In the running averaging method, the object computes an equally weighted average of a specified number of spectrum estimates defined by the SpectralAverages property. In the exponential method, the object computes the average over samples weighted by an exponentially decaying forgetting factor.

SpectralAverages — Number of spectral averages

8 (default) | positive integer

Number of spectral averages, specified as a positive integer. The spectrum estimator computes the current power spectrum or power density spectrum estimate by averaging the last N estimates. N is the number of spectral averages defined in the SpectralAverages property.

Dependencies

This property applies when you set AveragingMethod to 'Running'.

Data Types: double

ForgettingFactor — Forgetting factor

0.9 (default) | scalar in the range (0,1]

Specify the exponential weighting forgetting factor as a scalar value greater than zero and smaller than or equal to one.

Tunable: Yes

Dependencies

This property applies when you set AveragingMethod to 'Exponential'.

Data Types: single | double

ReferenceLoad — Reference load

1 (default) | positive scalar

Specify the load that the spectrum estimator uses as a reference to compute power values as a real, positive scalar in ohms.

Data Types: single | double

SidelobeAttenuation — Side lobe attenuation of window

60 dB (default) | positive scalar

Specify the side lobe attenuation of the window as a real, positive scalar, in decibels (dB).

Dependencies

This property applies when you set Method to 'Welch' and Window to 'Chebyshev' or 'Kaiser'.

Data Types: double

OutputMaxHoldSpectrum — Output max-hold spectrum

false (default) | true

Set this property to true so that the spectrum estimator computes and outputs the max-hold spectrum of each input channel. The max-hold spectrum is computed by keeping, at each frequency bin, the maximum value of all the power spectrum estimates.

OutputMinHoldSpectrum — Output min-hold spectrum

false (default) | true

Set this property to true so that the spectrum estimator computes and outputs the min-hold spectrum of each input channel. The min-hold spectrum is computed by keeping, at each frequency bin, the minimum value of all the power spectrum estimates.

SampleRate — Sample rate of input

1 (default) | scalar

Sample rate of the input in Hz, specified as a finite numeric scalar. The sample rate is the rate at which the signal is sampled in time.

Data Types: single | double

Usage

Syntax

pxx = SE(x)

[pxx,pmax] = SE(x)

[pxx,pmin] = SE(x)

[pxx,pmax,pmin] = SE(x)

Description

pxx = SE(x) computes the power spectrum or power-density spectrum, pxx, of the input signal, x. The System object treats the columns of x as independent channels.

[pxx,pmax] = SE(x) also computes the max-hold frequency spectrum, pmax, of x. To determine the max-hold spectrum, the method keeps the maximum of all the power spectrum estimates computed at each frequency bin. Set OutputMaxHoldSpectrum to true to obtain the max-hold spectrum.

[pxx,pmin] = SE(x) also computes the min-hold frequency spectrum, pmin, of x. To determine the min-hold spectrum, the method keeps the minimum of all the power spectrum estimates computed at each frequency bin. Set OutputMinHoldSpectrum to true to obtain the min-hold spectrum.

[pxx,pmax,pmin] = SE(x) computes the power spectrum or power-density spectrum, the max-hold spectrum, and the min-hold spectrum of x. Set OutputMaxHoldSpectrum and OutputMinHoldSpectrum to true to obtain the max-hold and min-hold spectra.

Input Arguments

x — Input signal

vector | matrix

Input signal, specified as a vector or matrix. The row length of x is the frame size or channel length. Each column of x is treated as a separate channel. The column length of x is the number of channels.

Data Types: single | double
Complex Number Support: Yes

Output Arguments

pxx — Power or power-density spectrum estimate

vector | matrix

Power or power-density spectrum estimate, returned as a vector or matrix of the same data type and complexity as the input signal, x.

When FFTLengthSource is set to:

• 'Auto' –– The size of pxx is same as the size of the input signal, x.

• 'Property' –– The size of pxx is the same as the specified FFT length.

By default, the unit of pxx is 'Watts'. You can also specify the spectrum to be in 'dBm' or 'dBW' through the PowerUnits property.

Data Types: single | double
Complex Number Support: Yes

pmax — Max-hold spectrum estimate

vector | matrix

Max-hold spectrum estimate, returned as a vector or matrix of the same size, data type, and complexity as the output signal, pxx.

Data Types: single | double
Complex Number Support: Yes

pmin — Min-hold spectrum estimate

vector | matrix

Min-hold spectrum estimate, returned as a vector or matrix of the same size, data type, and complexity as the output signal, pxx.

Data Types: single | double
Complex Number Support: Yes

Object Functions

To use an object function, specify the System object as the first input argument. For example, to release system resources of a System object named obj, use this syntax:

release(obj)

Specific to dsp.SpectrumEstimator

getFrequencyVector	Vector of frequencies at which estimation is done
getRBW	Resolution bandwidth of spectrum

Common to All System Objects

step	Run System object algorithm
release	Release resources and allow changes to System object property values and input characteristics
reset	Reset internal states of System object

Examples

Power Spectrum of Multichannel Sinusoidal Signal

Compute the power spectrum of a multichannel sinusoidal signal using the dsp.SpectrumEstimator System object™. You can get the vector of frequencies at which the spectrum is estimated using the getFrequencyVector function. To compute the resolution bandwidth of the estimate (RBW), use the getRBW function.

Generate a three-channel sinusoid sampled at 1 kHz. Specify sinusoidal frequencies of 100, 200, and 300 Hz. The second and third channels have their phases offset from the first by and , respectively.

sineSignal = dsp.SineWave(SamplesPerFrame=1000,SampleRate=1000, ...
    Frequency=[100 200 300],PhaseOffset=[0 pi/2 pi/4]);

Estimate and plot the one-sided spectrum of the signal. Use the dsp.SpectrumEstimator object for the computation and the dsp.ArrayPlot for the plotting.

estimator = dsp.SpectrumEstimator(FrequencyRange='onesided');
plotter = dsp.ArrayPlot(PlotType='Line',YLimits=[0 0.75], ...
    YLabel='Power Spectrum (watts)',XLabel='Frequency (Hz)');

Step through to obtain the data streams and display the spectra of the three channels.

y = sineSignal();
pxx = estimator(y);
plotter(pxx)

Get the vector of frequencies at which the spectrum is estimated in Hz, using the getFrequencyVector function.

f = getFrequencyVector(estimator);

Compute the resolution bandwidth (RBW) of the estimate using the getRBW function.

getRBW(estimator)

ans =

    0.0015

The resolution bandwidth of the signal power spectrum is 0.0015 Hz. This frequency is the smallest frequency that can be resolved on the spectrum.

Spectral Estimation Using Filter Bank

Compare spectral estimates of sinusoids embedded in white Gaussian noise using a Hann window-based Welch method and filter bank method.

Initialization

Initialize two dsp.SpectrumEstimator objects. Specify one estimator to use the Welch-based spectral estimation technique with a Hann window. Specify the other estimator to use an analysis filter bank to perform the spectral estimation. Specify a noisy sine wave input signal with 4 sinusoids at 0.16, 0.2, 0.205, and 0.25 cycles/sample. View the spectral estimate using an array plot.

FrameSize = 420;
Fs = 1;
sinegen = dsp.SineWave('SampleRate',Fs,...
    'SamplesPerFrame',FrameSize,...
    'Frequency',[0.16 0.2 0.205 0.25],...
    'Amplitude',[2e-5 1  0.05  0.5]);
NoiseVar = 1e-10;
numAvgs = 8;

hannEstimator = dsp.SpectrumEstimator('PowerUnits','dBm',...
    'Window','Hann','FrequencyRange','onesided',...
    'SpectralAverages',numAvgs,'SampleRate',Fs);

filterBankEstimator = dsp.SpectrumEstimator('PowerUnits','dBm',...
    'Method','Filter bank','FrequencyRange','onesided',...
    'SpectralAverages',numAvgs,'SampleRate',Fs);

spectrumPlotter = dsp.ArrayPlot(...
    'PlotType','Line','SampleIncrement',Fs/FrameSize,...
    'YLimits',[-250,50],'YLabel','dBm',...
    'ShowLegend',true,'ChannelNames',{'Hann window','Filter bank'});

Streaming

Stream the input. Compare the spectral estimates computed using the Hann window and the analysis filter bank

for i = 1:1000
    x = sum(sinegen(),2) + sqrt(NoiseVar)*randn(FrameSize,1);
    Pse_hann = hannEstimator(x);
    Pfb = filterBankEstimator(x);
    spectrumPlotter([Pse_hann,Pfb])
end

The Hann window misses the peak at 0.205 cycles/sample. In addition, the window has a significant spectral leakage that makes the peak at 0.16 cycles/sample hard to distinguish, and the noise floor is not correct.

The filter bank estimate has a very good resolution with no spectral leakage.

Power and Max-Hold Spectra of Noisy Sine Wave

Generate a sine wave.

sineWave = dsp.SineWave(Frequency=100,...
    SampleRate=1000, ...
    SamplesPerFrame=1000);

Use the spectrum estimator to compute the power spectrum and the max-hold spectrum of the sine wave. Use the Array Plot to display the spectra.

SE = dsp.SpectrumEstimator(...
    SampleRate=sineWave.SampleRate,...
    SpectrumType='Power',PowerUnits='dBm', ...
    FrequencyRange='centered',...
    OutputMaxHoldSpectrum=true);
plotter = dsp.ArrayPlot(PlotType='Line',...
    XOffset=-500, ...
    YLimits=[-60 30],...
    Title='Power Spectrum of 100 Hz Sine Wave', ...
    YLabel='Power Spectrum (dBm)',...
    XLabel='Frequency (Hz)');

Add random noise to the sine wave. Stream in the data, and plot the power spectrum of the signal.

for ii = 1:10
    x = sineWave() + 0.05*randn(1000,1);
    [Pxx,Pmax] = SE(x);
    plotter([Pxx Pmax])
end

Algorithms

Welch's Method of Averaged Modified Periodograms

When you choose the Welch method, the power spectrum estimate is averaged modified periodograms.

Given the signal input x:

1. Multiply x by the window and scale the result by the window power.

2. Compute the FFT of the signal Y and take the square magnitude using Z = Y.*conj(Y).

3. Compute the current power spectrum estimate by taking the moving average of the last N number of Zs, and scaling the answer by the sample rate. For more information on the moving average methods, see Averaging Method.

Filter Bank

When you select the filter bank method, the spectrum estimator uses an analysis filter bank to estimate the power spectrum. The filter bank splits a broadband input signal, x(n), of sample rate fs into multiple narrow band signals, y₀(m), y₁(m), … , y_M-1(m), of sample rate fs/M.

The variable M represents the number of frequency bands in the filter bank. When you specify FFT length, M equals the FFT length. When you do not specify FFT length, M is equal to the number of rows in the input signal. The number of taps per frequency band sets the number of filter coefficients for each frequency band of the filter bank. The total number of filter coefficients is equal to the number of taps per band multiplied by the number of frequency bands M. For more information on the analysis filter bank and how it is implemented, see the More About and the Algorithm sections in dsp.Channelizer.

After the broadband input signal is split into multiple narrow bands, the spectrum estimator computes the power in each narrow band using the following equation. Each Z_i value becomes the estimate of the power over that narrow frequency band.

$$Z_i=\frac{1}{L}{\displaystyle \sum _{m=0}^{L-1}{\left|y_i[m]\right|}^2}$$

where L is the length of the narrow band signal y_i(m) and i = 1, 2, …, M−1.

The power values in all the narrow bands (denoted by Z_i) form the Z vector.

$$Z=[Z_0,Z_1,Z_2,\cdots ,Z_{M-1}]$$

The filter bank estimator algorithm averages the current Z vector with the previous Z vectors using one of the two moving average methods: running or exponential weighting. The output of the averaging operation forms the spectral estimate vector. For more information on the two averaging methods, see Averaging Method.

Averaging Method

The algorithm calculates the moving average using one of these methods:

• Running — For each frame of input, average the last N-scaled Z vectors, which are computed by the algorithm. The variable N is the value you specify for the number of spectral averages. If the algorithm does not have enough Z vectors, the algorithm uses zeros to fill the empty elements.

  

• Exponential — The moving average algorithm using the exponential weighting method updates the weights and computes the moving average recursively for each Z vector that comes in by using the following recursive equations:

  $$\begin{array}{l}{w}_N=\lambda w_{N-1}+1\\ {\overline{z}}_N=\left(1-\frac{1}{w_N}\right){\overline{z}}_{N-1}+\left(\frac{1}{w_N}\right)z_N\end{array}$$

  • λ — Forgetting factor

  • $$w_N$$ — Weighting factor applied to the current Z vector

  • $$z_N$$ — Current Z vector

  • $${\overline{z}}_{N-1}$$ — Moving average until the previous Z vector

  • $$\left(1-\frac{1}{w_N}\right){\overline{z}}_{N-1}$$ — Effect of the previous Z vectors on the average

  • $${\overline{z}}_N$$ — Moving average including the current Z vector

References

[1] Hayes, Monson H. Statistical Digital Signal Processing and Modeling. Hoboken, NJ: John Wiley & Sons, 1996

[2] Kay, Steven M. Modern Spectral Estimation: Theory and Application. Englewood Cliffs, NJ: Prentice Hall, 1999

[3] Stoica, Petre and Randolph L. Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005

[4] Welch, P. D. “The use of fast Fourier transforms for the estimation of power spectra: A method based on time averaging over short modified periodograms,” IEEE Transactions on Audio and Electroacoustics, Vol. 15, 1967, pp. 70–73.

[5] Harris, F.J. Multirate Signal Processing for Communication Systems. Prentice Hall. 2004, pp. 208–209.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

Usage notes and limitations:

• System Objects in MATLAB Code Generation (MATLAB Coder)

Version History

Introduced in R2013b