Simple but Effective Spectrum Averaging

Digitally modulated signals and Gaussian noise are random signals, and randomness is apparent in their spectra.  For example, Figure 1 (left) shows the power spectral density (PSD) of a QAM-16 signal with added Gaussian noise.  While this PSD is perfectly correct, it has a large variation around the average values of the signal and the noise, making it difficult to estimate levels.  We can reduce the variation by averaging the power density of several PSD’s [1,2].  An averaged PSD is shown in Figure 1 (right).

In PSD averaging, we compute several independent PSD’s.  Suppose a given frequency bin k of the first PSD has power P₁(k), the same bin of the 2^nd PSD has power P₂(k), etc.  Given M PSD’s, averaging P₁(k), P₂(k), … P_M(k) gives the average power of bin k.  Averaging all the bins in this manner gives us the average PSD.  Thus, each point in the right plot of Figure 1 is an average.  Keep in mind that this is not an average over frequency; rather, it is an average of each frequency bin over multiple PSD’s.

This article is available in PDF format for easy printing.

In this article, I provide a Matlab function that performs exponential PSD averaging, using first-order infinite impulse response (IIR) filtering [3] to continuously average the PSD bins.  This approach works well for computing the spectrum of a long-duration signal over time, because the spectrum is constantly updated as new PSD’s are computed.  Conveniently, the time constant of the PSD averaging is determined by the single adjustable parameter α.  I also provide a Matlab function for conventional (unweighted) PSD averaging.  Neither function requires any canned code other than the Fast Fourier Transform (FFT), although I do use the Matlab hann window function for convenience.

Figure 2 is a simplified diagram of the averaging process.  First, we acquire N samples of x(n).  We then compute the N-point Discrete Fourier Transform (DFT).  The DFT is shown as X(k), and the PSD is proportional to |X(k)|² [4].  A new PSD is computed every N samples of x.  We filter each PSD bin over all the PSD’s, which gives us the averaged PSD.

Figure 1.  Left:  Power Spectral Density (PSD) of QAM-16 signal with noise.

            Right:  Average of multiple PSD’s.

Figure 2. Simplified Block Diagram of Spectrum Averaging.

Matlab Function psd_exp

The function psd_exp performs exponential PSD averaging, using first-order IIR filters to average the PSD frequency bins.  A listing is provided in Appendix A of the PDF version.  The function call is as follows:

    [PdB,f] = psd_exp(x,N,alpha,fs)

where

x = input signal vector

N = length of DFT, samples

alpha = averaging parameter, 0 < alpha <= 1

fs = sample frequency, Hz

PdB = averaged power spectral density vector of length N/2, dBW/bin

f = frequency vector of PdB, 0 <= f < fs/2, Hz

The spectra of contiguous length-N segments of x are averaged.  Letting L be the length of x, L must be greater than or equal to N.  L is typically an integer multiple of N.  For example, if L = 32*N, then 32 spectra are averaged.  

Alpha must be greater than zero and less than or equal to one.  For alpha = 1, no spectrum averaging occurs.  1/alpha is approximately equal to the time constant of the exponential averagers.  For example, setting alpha = 1/8 results in an averaging time constant of roughly 8 PSD’s.  Lower alpha provides more smoothing of the spectrum.  To allow the averaged PSD to settle, L/N = number of PSD’s should normally be greater than 2/alpha.

The PSD output PdB has units of dBW/bin, assuming a 1-ohm load.  Besides outputting the final averaged PdB, the function plots each PdB vector as it is calculated, showing the effect of averaging over time.

Example:  Sine Wave Plus Gaussian Noise

Here we compute the PSD of a sinewave plus Gaussian noise using psd_exp.  Sine frequency is 9 Hz and sample rate is 100 Hz.  We set the sine amplitude A = sqrt(2) volts, which gives power of A²/2 = 1 watt into 1 ohm.  The following Matlab code generates the sinewave and adds Gaussian noise:

    fs = 100;             % Hz  sample frequency
    Ts = 1/fs;            % s  sample time
    L = 2^14;             % samples  length of signal x
    f0 = 9;               % Hz  sine frequency
    A = sqrt(2);          % sine amplitude
    n= 0:L-1;             % samples  time index
    u = A*sin(2*pi*f0*n*Ts);    % sinewave
    x = u + 0.016*randn(1,L);    % sine + Gaussian noise

Next, we compute the non-averaged PSD by setting alpha = 1.  DFT length is 512.  In the following code, we use only the first 512 samples of the signal x, so only one DFT is performed.

    N = 512;                  % samples  length of DFT
    alpha = 1;                % no averaging case
    [PdB1,f]= psd_exp(x(1:N),N,alpha,fs);    % dBW/bin  PSD

The resulting spectrum of N/2 frequency samples is shown in Figure 3 (left), where we see a large variation in the noise floor.  Note that the sine amplitude is 1.8 dB less than 0 dBW (1 watt).  This error is due to the processing loss [5] of the Hanning window function used by psd_exp.  The error only occurs for narrow spectral components that occupy at most a few bins.  For wider spectra, such as that of Figure 1, the displayed dBW/bin is correct.

Finally, we find the averaged spectrum using all 2^14 samples of x, with alpha = 1/8.  This alpha corresponds to an averaging time constant of roughly 1/alpha = 8.  Given the DFT length of 512, we average L/512 = 2^14/512 = 32 DFT’s.

    alpha = 1/8;
    [PdB2,f]= psd_exp(x,N,alpha,fs);    % dBW/bin  average PSD

The spectrum is shown in Figure 3 (right).  The variation of the noise floor has been reduced by spectrum averaging.  Again, this result was achieved by averaging each frequency bin over all the PSD’s.  From the plot, we can estimate the average noise spectral density to be about -60 dBW/bin.

It is instructive to plot the PSD of a single bin of the PSD vs. the sequence number i of the PSD, with and without averaging.  Figure 4 shows this for a bin in the noise floor with bin index k = 120 (f = k*fs/N = 25.39 Hz).  It requires at least 2/alpha = 16 PSD’s for the averaged PSD (red curve) to settle.  Compared to the unaveraged case (blue curve), we see that the PSD variation of this single bin is reduced significantly.

Figure 3.  PSD of sinewave plus Gaussian noise.

Left:  no averaging.  Right:  with averaging, alpha = 1/8.

Figure 4.  PSD of a single bin, with and without averaging.

Blue:  no averaging.  Red:  with averaging, alpha = 1/8.

i = sequence number of PSD.  bin index k = 120 (f = 25.39 Hz).

How it Works

Figure 5 shows the sequence of operations of Matlab function psd_exp, whose code is listed in Appendix A of the PDF version.  The code is not long, but it does rely on several DSP concepts:  DFT, windowing for the DFT, PSD, and exponential averaging/first-order IIR filters.  Most of the formulas shown below were developed in earlier posts.  I use lower-case names for time signals and upper-case names for frequency spectra.  Thus “x_win” is a time signal and its DFT “X” is a frequency spectrum.  The operations are described as follows.

Acquisition:  Acquire N samples of the input signal x(n).This signal vector is called xi.  In the code, the captured segments are contiguous, although this may not be required to obtain a valid spectrum.

Windowing:  Windowing multiplies each sample of xi by each element of a window function of the same length.  The window function used is the simple Hanning, or Hann window [6]: 

$$ win(n)= 0.5(1 - cos(2\pi n/N)), \quad n = 0:N-1 \qquad (1)$$

In the code, Equation 1 is implemented by the Matlab function hann(N,’periodic’)’.  (The apostrophe at the end of the function call converts the window’s column vector to a row vector).  If lower spectral sidelobes are needed, a Kaiser window function can be used [7 - 9].

The amplitude of the window function will affect the spectrum amplitude if not compensated.  The scaling factor [4] to normalize the window's average power to 1 watt is: 

$$ scaling \,factor = \sqrt{ \frac{N}{\sum_{n=0}^{N-1}(win(n))^2} } \qquad (2) $$

If absolute power values are not required, the scaling can be omitted.

DFT:  The N-point FFT of x_win is computed, and the first N/2 samples are retained in the one-sided spectrum X.

PSD:  The PSD [4] is calculated using:

$$ PX = \frac{2}{N^2}|X(k)|^2, \quad k = 0:N/2 -1 \qquad watts/bin \qquad (3) $$

where k is frequency index (frequency bin) and the vector PX has units of watts/bin, assuming a 1-ohm load.  For power into a load of R ohms, just divide the result by R.  If absolute power values are not required, the scaling by 2/N² can be omitted.

Exponential Averaging:  The difference equation of an exponential averager or first-order lowpass IIR filter [3] is:

$$ y(n) = \alpha*x(n) + (1-\alpha)*y(n-1) \qquad (4) $$

where input x(n) is a general time signal (not related to x in the code) and 0 < α < 1.  To average our PSD PX, each bin k of PX is applied to an exponential averager.  So, using explicit subscripts:

$$ PY(k,i) = \alpha*PX(k,i) + (1-\alpha)*PY(k,i-1), \quad k= 0:N/2-1 \qquad (5) $$

where PY is the averaged PSD and i is the sequence number of PX and PY (index of the FOR loop).  Equation 5 is evaluated N/2 times for each i.  In the Matlab code, PX and PY are vectors in k of length N/2 and the subscripts k and i are implicit.  Note that for α = 1, PY(k,i) = PX(k,i) and no averaging occurs.

Convert to dB:  PY is converted from watts/bin to dBW/bin using:

$$ PdB = 10*log_{10}(PY) \quad dBW/bin \qquad (6)  $$

Where PdB is a vector of length N/2.

Plotting the PSD:  The function plots each PdB vector as it is calculated, showing the effect of averaging over time.  (There is a brief pause in the code after each plot).  The dB-axis of the plot is fixed, so the input signal x may require scaling to plot properly (or the axis scaling can be changed in the code).

Conventional PSD Averaging

A Matlab function for conventional PSD averaging is listed in Appendix B of the PDF version.  It uses all the above operations, except it substitutes unweighted averaging for exponential averaging.  Let the signal x(n) have length L and let the DFT length be N.  Then the number of PSD’s averaged is:

$$ Navg = fix(L/N) \qquad (7) $$

To average our PSD PX, each bin k of PX is summed over all PSD’s.  Using explicit subscripts, the running sum of the PSD’s is:

$$ Psum(k,i)= PX(k,i) + Psum(k,i-1), \quad k= 0:N/2-1 \qquad (8) $$

where i is the sequence number of PX and Psum (Index of the FOR loop).  The average power of the PSD is then:

$$ PY(k) = \frac{Psum(k,Navg)}{Navg}, \quad k = 0:N/2-1 \qquad (9) $$

In the Matlab code (Appendix B of the PDF version), PX, Psum, and PY are vectors in k of length N/2 and the subscripts are implicit.

Figure 5.  Sequence of operations of Matlab function psd_exp.

---

For Appendices, see the PDF version.

---

References

1.  Lyons, Richard G.Understanding Digital Signal Processing, 3^rd Ed., Pearson, 2011, Sections 11.2 and 11.6.3.

2.  Robertson, Neil, “Use Matlab Function pwelch to Find Power Spectral Density – or Do It Yourself,” DSPRelated.com, Jan., 2019, https://www.dsprelated.com/showarticle/1221.php

3.  Robertson, Neil, “The First-Order IIR Filter – More Than Meets the Eye,” DSPRelated.com, Nov., 2025,

https://www.dsprelated.com/showarticle/1769.php

4.  Robertson, Neil, “The Power Spectrum,” DSPRelated.com, Oct., 2016, https://www.dsprelated.com/showarticle/1004.php

5.  Robertson, Neil, “Evaluate Window Functions for the Discrete Fourier Transform,” DSPRelated.com, Dec., 2018, https://www.dsprelated.com/showarticle/1211.php

6.  Lyons, p 91.

7.  The Mathworks website, “kaiser”, https://www.mathworks.com/help/signal/ref/kaiser.html

8.  Robertson, Neil, “A Simplified Matlab Function for Power Spectral Density,” DSPRelated.com, March, 2020, https://www.dsprelated.com/showarticle/1333.php

9.  Lyons, p 197 – 201.

                               Neil Robertson     January 2026

Sign in