easyspec: A spectrum-analyzer like plotter for MATLAB

I’m using MATLAB since about 1992. One of the things that I’ve always missed, as a communication engineer, is a quick and dirty spectrum plot of signals. I want to see what the spectrum analyzer will show me when I’ll run the signal through it, and not some analytic plot, which looks nothing like.

The “easyspec” MATLAB function does exactly that. Using this function on a time signal has the feel of plugging in a spectrum analyzer: You get some simple, concise image (scaled in dBs, for us engineers), which is as reliable as if you measured the signal. Meaning, that if the signal is not stationary (choose whatever sense you want for “stationary”) you may want to re-run the function to see if you get the same result.

Let’s try a simple session (you’ll need to copy the “easyspec” function from the bottom of this page and put it in your working directory or a directory in the MATLAB search path):

>> t=0:1/44100:10;
>> s=0.25*sin(2*pi*440*t) + sin(2*pi*1000*t) + randn(1,441001);
>> easyspec(s, 44100);

This shows the spectrum of two tones, one at 440Hz, and another at 1000Hz, sampled at 44.1 kHz. Note that since we hand over the sample rate to easyspec, it will show the real-life frequency on the x axis.

The special thing about this function, is that it chooses the FFT segments at random. As a result, subsequent calls to this function on the same signal will result in slightly differently plots. This is exactly what we expect from a real-life spectrum analyzer.

As given below, this function picks (at random) 100 segments of 4096 samples, applies a Hanning window to each, and pads the vector to a length of 16384. The absolute value of the FFT of this vector is calculated.

What is displayed is the average of these 100 absolute value vectors, normalized so that a sine wave of amplitude 1 will show as approximately 0 dB on the spectrum plot. This is in the spirit of the “video average” function of common spectrum analyzers, but it works differently (spectrum analyzers usually average on the dB values). Anyhow, this makes a smoother display of truly random signals (such as noise).

The most common question I get is where the 17.127 constant came from. The honest answer is that I simply played with the number until I got 0 dB for a unity sine wave. It’s true that I could have made some calculations to reach a constant with theory behind, but I found it pretty pointless, since the accuracy is compromised anyhow. (The main reason is that the reading depends on the frequency, as each FFT bin’s response to a pure sine wave is the Hanning window’s Fourier transform).

The reason for picking the segments at random comes from real-life spectrum analyzers: Almost always, there is no synchronization between the sweep instances and the measured signal. Also, let’s keep in mind that the term “spectrum” comes from the world of probability theory, and it’s the Fourier Transform of the autocorrelation function of a side-sense stationary random signal. Now we may ask ourselves what meaning there is to a spectrum of a deterministic signal. What’s random about a plain sine wave that we apply to a spectrum analyzer?

The answer is that the time segment, on which we measure the signal, is random. If it’s not, as in MATLAB’s native “spectrum” function, a misleading result may be displayed as a result of some relation between the measured signal’s periodicity and the jumps between the FFT segments. This doesn’t happen when the signal is simple noise, but as soon as there is anything periodic about it, strange things can happen.

And finally, here is the function in question. Copy and save it as ‘easyspec.m’. It’s released to the public domain:

function [s,f]=easyspec(x,fs)
%EASYSPEC Easy plot of spectrum estimate
%         S=EASYSPEC(X) will return a spectrum vector of
%         X. [S,F]=EASYSPEC(X,FS) will also return a freuqency
%         axis in F, while the sample frequency is given in FS.
%         EASYSPEC(X) and EASYSPEC(X,FS), will plot the spectrum
%         and return the S vector in MATLAB's "ans".
%         Notes:
%         1. Spectrum values are in dB.
%         2. The frequency domain goes from DC to the Nyquist 
%            frequency. The "mirror part" of the spectrum is
%            omitted.
%         3. The sample segments, from which the spectrum is
%            calculated are picked by random. The function might
%            return significantly different results each time it's
%            called to if X isn't stationary.
%         4. EASYSPEC uses a hanning window and zero-pads by a
%            factor of 4. The spectrum vector will look smooth.

% Eli Billauer, 17.1.01
% This function is released to the public domain; Any use is allowed.

if nargin==0
  error('No input vector given');
end

if (nargin==1)
  fs=2;
end

NFFT=16384; NWIN=NFFT/4;
LOOP=100;
win=hanning(NWIN)';

x=x(:)'*(17.127/NFFT/sqrt(LOOP));
n=length(x);
maxshift=n-NWIN;

if (n<2*NWIN)
  error(['Input vector should be at least of length '...
      num2str(2*NWIN)]);
end

s=zeros(1,NFFT);

for i=1:LOOP
  zuz=floor(rand*maxshift);
  s=s+abs(fft([win.*x(1+zuz:NWIN+zuz) zeros(1,NFFT-NWIN)])).^2;
end

s=10*log10(s(1:NFFT/2));
f=linspace(0,fs/2,NFFT/2);

if nargout==0
  hold off;
  plot(f,s);
  ylabel('Power Spectrum [dB]');
  xlabel('Frequency');
	grid on; zoom on;
end