I am trying to reproduce a numerical model of classical spectrum analyzer, described in "Fundamentals of Spectrum Analysis" by Christoph Rauscher.

Classical model consists of a few local oscillators (by number of frequency sweep points). The LO signals should be mixed with analyzed signal one-by-one, and then products should be convoluted with low pass filter impulse response. In the end, power of each filtered signal should be calculated, thus one can plot these values vs. swept frequencies.

In spectrum analysis the most popular IF band pass filter is Gaussian, so i use Gaussian low pass as equivalent.

Power spectrum in dB, normalized, 10 MHz sampling, 1 MHz carrier, 100 kHz span, 501 sweep points, 10 kHz RBW, 1 ms sweep time

The resulting spectrum looks bad (see figure). 3 dB bandwidth is not reproduced. Gaussian form distorted near carrier frequency.

Any suggestions?

Matlab MWE:

clear variables;
close all;

%% Signal properties
fc = 1e6;                               % Carrier frequency, Hz
fs = 10e6;                              % Sampling frequency, Hz

%% Spectrum analyzer specification
span = 100e3;                           % Span window, Hz
center = 1e6;                           % Center frequency, Hz
sweep = 1e-3;                           % Sweep time, s
rbw = 10e3;                             % Resolution filter bandwidth, Hz
points = 501;                           % Sweep points

%% Signal synthesis
samples = sweep*fs+1;

% Time scale
time = linspace(0, sweep, samples);

% The signal
signal = cos(2*pi*fc*time);

%% Gaussian filter design
order = ceil(fs/2/rbw)*10;
alpha = 10;
cutoff = rbw/2;
window = gausswin(order+1, alpha);
b = fir1(order, cutoff/(fs/2), 'low', window, 'scale');

%% Signal analysis
start = center - span/2;
stop = center + span/2;
freqs = linspace(start, stop, points);
power = 10*log10(sum(signal.^2)/samples*(samples-1-order))-3;

product = zeros(samples, points);
values = zeros(1,points);
for i = 1:points
    product(:,i) = (signal.*cos(2*pi*freqs(i)*time))';
    product(:,i) = filter(b,1,product(:,i));
    values(i) = 10*log10(sum(product((order+1):end,i).^2))-power;

xlim([start stop]);
ylim([-120 0]);
grid on;

1 Answer 1


I believe the issue is the OP is simulating a spectrum analyzer with a "Zero-IF" architecture so is seeing the image issues as the sweep signal approach the input signal (since both are real signals). The product is down-converted to baseband with a real LO and passed through a low-pass Gaussian filter, so the resulting signal as it approach the DC center as a very low frequency real sine-wave has increasingly fewer cycles (and ultimately partial cycles) which then contributes to high variability in the power estimate. By using a complex LO instead, the envelope of the resulting signal is constant eliminating that variability (or by observing the two real components of the baseband complex signal $e^{j\omega t} = \cos(\omega t) + j\sin(\omega t)$ we see that as one partial cycle is decreasing the other is increasing to maintain the constant power at the detector).

I confirmed this was the case by changing the LO to a complex signal in the OP's code which would then match the traditional architecture of an actual Zero-IF receiver:

product(:,i) = (signal.*exp(-j*2*pi*freqs(i)*time))';

and use a complex conjugate product for the power (which also could be done using abs(x).^2)

values(i) = 10*log10(sum(product((order+1):end,i).*(conj(product((order+1):end,i)))))-power;

Resulting in:

clean result

The other option (that I would NOT recommend for simulation! - unless you need to prove to yourself that it is identical) is to move the filter from being a low-pass filter to a band-pass filter and offset the LO such that an IF frequency passes through the filter (matching traditional spectrum analyzer architectures).

For simulation purposes I would always simulate the equivalent baseband signal using complex signals rather than actual carriers which would be significantly more time consuming to process.

Further, this process can be done in one operation using an FFT which would be significantly faster and for a stationary input such as this produce an identical result. The Gaussian filter could be implemented by multiplying the time domain waveform by the Gaussian window prior to taking the FFT. For non-stationary signals further additions could be made to properly emulate the effect of the sweep time.

As far as the bandwidth, the Gaussian window used in the OP's code with $\alpha = 10$ has an equivalent noise BW of 5.64 bins as computed from:

$$\text{ENBW} = N\frac{\sum (w[n]^2)}{(\sum w[n])^2} \tag{1} \label{1}$$

That said for a sampling rate of 10e6, and a target resolution BW of 10KHz, this would require a time capture of

$$N = 5.64 \times (10e6/10e3) = 5640 \text{ samples}$$

This would then also match the duration for the Gaussian filter for a 10KHz resolution BW in the OP's code by selecting the coefficients directly from the window function of that length, as in:

b =  gaussian(5640, 10)

I use equivalent noise bandwidth (ENBW) as resolution BW instead of the 3 dB BW, as it can be directly computed from the window as in \ref{1} and accurately indicates the same power that would be in a brick-wall filter of that bandwidth for a white noise source. The ENBW is slightly higher than the 3 dB bandwidth (on the order of 5% but depends on the window).

As recently reiterated by fred harris at the DSP Online Conference (https://www.dsponlineconference.com/), the Gaussian window is NOT the best choice. As he detailed, the reason Gaussian is often selected for spectral analysis is that it is believed to have the minimum time-bandwidth product, so minimum localization in time and frequency. However this is only true for a true-Gaussian which requires an infinite time support and only in that case when $t \rightarrow \pm\infty$ results in the minimum achievable time-bandwidth product BT of $BT = 0.5$. When we truncate the Gaussian response for the time-limited window (as required for practical implementation) the truncated Gaussian $BT>0.5$ and is no longer the minimum BT window function available for the case of finite time duration.

A window that has minimum time BW product under finite time support is the DPSS (Digital Prolate-Spheroid Sequence) window and very close to that is the Kaiser Window, both of which are available in MATLAB, Octave and Python scipy.signal. These windows also have decreasing stop-band attenuation which is also attractive for many applications. For more info on this see: https://ccrma.stanford.edu/~jos/sasp/Slepian_DPSS_Window.html , https://www.dsprelated.com/freebooks/sasp/Kaiser_DPSS_Windows_Compared.html#fig:dpsstest and https://ieeexplore.ieee.org/document/543677

  • $\begingroup$ Thank you very much! $\endgroup$
    – akj
    Sep 30, 2020 at 5:53
  • $\begingroup$ I thought that it works just as you say. The energy in a few samples (compared to very low frequency of fundamental component) in product is just not enough to complete PSD-curve. Also i thought that Gaussian window is used due to its invariant under the Fourier transform. $\endgroup$
    – akj
    Sep 30, 2020 at 6:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.