5.6.3. How a spectrum is computed in PyRPL¶

5.6.3.1. Definitions¶

We can show that the Fourier transform of the product is the convolution of the Fourier Transforms, such that:

Y[r] = 1/N sum_r’ X[r’] W[r-r’] (1)

To make sure the windowing function is well normalized, and to define the noise equivalent bandwidth of a given windowing function, we will study the 2 limiting cases where the initial time series is either a sinusoid or a gaussian distributed white noise.

5.6.3.2. Sinusoidal input¶

To simplify the calculations, we assume the period of the sinusoid is a multiple of the sampling rate:

x[k] = cos[2 pi m k/N]

= 1/2 (exp[i 2 pi m k/N] + exp[-2 pi i (N - m) k/N])

We obtain the Fourier transform:

X[r] = N/2 (delta[r-m] + delta[r-(N-m)]).

We deduce using (1), that the estimated spectrum is:

Y[r] = 1/2 (W[r - m] + W[r - (N-m)])

With the discrete fourier transform convention used here, we need to pay attention that the DC-component is for r=0, and the ?negative frequencies? are actually located in the second half of the interval [N/2, N]. If we take the single sided convention where the negative frequency side is simply ignored, the correct normalization in terms of V_pk (for which the maximum of the spectrum corresponds to the amplitude of the sinusoid) is the one for where max(W[r]) = 2.

Moreover, a reasonable windowing function will only have non-zero Fourier components on the few bins around DC, such that if we measure a pure sinusoid with a frequency far from 0, there wont be any significant overlap between the two terms, and we will measure 2 distinct peaks in the positive and negative frequency regions, each of them with the shape of the Fourier transform of the windowing function. Since the maximum of W[r] is located in r=0, we finally have:

sum_k w[k] = 2

5.6.3.3. White noise input¶

Once the normalization of the filter window has been imposed by the previous condition, we need to define the bandwidth of the window such that noise measurements integrated over frequency give the right variances.

Let’s take a white noise of variance 1.

<x[k] x[k’]> = delta(k-k’).

We would like the total spectrum in units of Vrms^2/Hz, integrated from 0 to Nyquist frequency to yield the same variance of 1. This is ensured by the Equivalent noise bandwidth of the filter window. To convert from V_pk^2 to V_rms^2/Hz, the spectrum is divided by the residual bandwidth of the filter window.

Let’s calculate:

sum_r <|Y[r]|^2> = (...) = N sum_k w[k]^2 <|x[k]|^2>

If we remind that x[k] is a white noise following <|x[k]|^2> = 1, we get:

sum_r <|Y[r]|^2> = N sum_k w[k]^2

So, since we want:

sum_r <|Z[r]|^2> df = 2, (indeed, we want to work with single-sided spectra, such that integrating over positive frequencies is enough)

with df the frequency step in the FFT, we need to choose:

rbw = N sum_k w[k]^2 df /4

In order to use dimensionless parameters for the filter windows, we can introduce the equivalent noise bandwidth:

ENBW = sum_k w[k]^2/(sum_k w[k])^2 = 1/4 sum_k w[k]^2

Finally, we get the expression of the rbw:

rbw = sample_rate ENBW

5.6.3.4. IQ mode¶

In iq mode, the signal to measure is fed inside an iq module, and thus, multiplied by two sinusoids in quadrature with each other, and at the frequency center_freq. The resulting I and Q signals are then filtered by 4 first order filters in series with each other, with cutoff frequencies given by span. Finally, these signals are measured simultaneously with the 2 channels of the scope, and we form the complex time serie c_n = I_n + i Q_n. The procedure described above is applied to extract the periodogram from the complex time-serie.

Since the data are complex, there are as many independent values in the FFT than in the initial data (in other words, negative frequencies are not redundant with positive frequency). In fact, the result is an estimation of the spectrum in the interval [center_freq - span/2, center_freq + span/2].

5.6.3.5. Baseband¶

In baseband mode, the signal to measure is directly fed to the scope and the procedure described above is applied directly. There are 2 consequences of the fact that the data are real:

1. The negative frequency components are complex conjugated (and thus redundant) wrt the positive ones. We thus throw away the negative frequencies, and only get a measurement on the interval [0, span/2]
2. The second scope channel can be used to measure another signal.

It is very interesting to measure simultaneously 2 signals, because we can look for correlations between them. In the frequency domains, these correlations are most easily represented by the cross-spectrum. We estimate the cross-spectrum by performing the product conjugate(fft1)*fft2, where fft1 and fft2 are the DFTs of the individual scope channels before taking their modulus square.

Hence, in baseband mode, the method curve() returns a 4x2^13 array with the following content: - spectrum1 - spectrum2 - real part of cross spectrum - imaginary part of cross spectrum

5.6.3.6. Proposal for a cleaner interface for spectrum analyzer:¶

To avoid baseband/2-channels acquisition from becoming a big mess, I suggest the following:

• The return type of the method curve should depend as little as possible from the particular settings of the instrument (channel2_baseband_active, display_units). That was the idea with scope, and I think that makes things much cleaner. Unfortunately, for baseband, making 2 parallel piplines such as curve_iq, curve_baseband is not so trivial, because curve() is already part of the AcquisitionModule. So I think we will have to live with the fact that curve() returns 2 different kinds of data in baseband and iq-mode.
• Moreover, in baseband, we clearly want both individual spectra + cross-spectrum to be calculated from the beginning, since once the abs() of the ffts is taken, it is already too late to compute conjugate(fft1)*fft2
• Finally, I suggest to return all spectra with only one “internal unit” which would be V_pk^2: indeed, contrary to rms-values unittesting doesn’t require any conversion with peak values, moreover, averaging is straightforward with a quadratic unit, finally, .../Hz requires a conversion-factor involving the bandwidth for unittesting with coherent signals

I suggest the following return values for curve():

• In normal (iq-mode): curve() returns a real valued 1D-array with the normal spectrum in V_pk^2
• In baseband: curve() returns a 4xN/2-real valued array with (spectrum1, spectrum2, cross_spectrum_real, cross_spectrum_imag). Otherwise, manipulating a complex array for the 2 real spectra is painful and inefficient.

Leo: Seems okay to me. One can always add functions like spectrum1() or cross_spectrum_complex() which will take at most two lines. Same for the units, I won’t insist on rms, its just a matter of multiplying sqrt(1/2). However, I suggest that we then have 3-4 buttons in the gui to select which spectra and cross-spectra are displayed.

Yes, I am actually working on the gui right now: There will be a baseband-area, where one can choose display_input1_baseband, input1_baseband, display_input2_baseband, input2_baseband, display_cross_spectrum, ‘display_cross_spectrum_phase’. And a “iq-area” where one can choose center_frequency and input. I guess this is no problem if we have the 3 distinct attributes input, input1_baseband and input2_baseband, it makes thing more symmetric...

5.6.3.7. IQ mode with proper anti-aliasing filter¶

When the IQ mode is used, a part of the broadband spectrum of the two quadratures is to be sampled at a significantly reduced sampling rate in order to increase the number of points in the spectrum, and thereby resolution bandwidth. Aliasing occurs if significant signals above the scope sampling rate are thereby under-sampled by the scope, and results in ghost peaks in the spectrum. The ordinary way to get rid of this effect is to use excessive digital low-pass filtering with cutoff frequencies slightly below the scope sampling rate, such that any peaks outside the band of interest will be rounded off to zero. The following code implements the design of such a low-pass filter (we choose an elliptical filter for maximum steepness):

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# the overall decimation value
decimation = 8

# elliptical filter runs at ell_factor times the decimated scope sampling rate
ell_factor = 4

wp = 0.8/ell_factor # passband ends at xx% of nyquist frequency
ws = 1.0/ell_factor # stopband starts at yy% of nyquist frequency
gpass = 5. # jitter in passband (dB)
gstop = 20.*np.log10(2**14)  # attenuation in stopband (dB)
#gstop = 60  #60 dB attenuation would only require a 6th order filter
N, Wn = signal.ellipord(wp=wp, ws=ws, gpass=gpass, gstop=gstop, analog=False)  # get filter order
z, p, k = signal.ellip(N, gpass, gstop, Wn, 'low', False, output='zpk')  # get coefficients for implementation
b, a = signal.ellip(N, gpass, gstop, Wn, 'low', False, output='ba')  # get coefficients for plotting
w, h = signal.freqz(b, a, worN=2**16)
ww = np.pi / 62.5  # scale factor for frequency axis (original frequency axis goes up to 2 pi)

# extent w to see what happens at higher frequencies
w = np.linspace(0, np.pi, decimation/ell_factor*2**16, endpoint=False)
# fold the response of the elliptical filter
hext = []
for i in range(decimation/ell_factor):
    if i%2 ==0:
        hext += list(h)
    else:
        hext += reversed(list(h))
h = np.array(hext)
# elliptical filter
h_abs = 20 * np.log10(abs(h))

# 4th order lowpass filter after IQ block with cutoff of decimated scope sampling rate
cutoff = np.pi/decimation
butter = 1.0/(1.+1j*w/cutoff)**4
butter_abs = 20 * np.log10(abs(butter))

# moving average decimation filter
M = float(decimation) # moving average filter length
mavg = np.sin(w*float(M)/2.0)/(sin(w/2.0)*float(M))
mavg_abs = 20 * np.log10(abs(mavg))

# plot everything together and individual parts
h_tot = h_abs + mavg_abs + butter_abs
plt.plot(w/ww, h_tot, label="all")
plt.plot(w/ww, h_abs, label="elliptic filter")
plt.plot(w/ww, butter_abs, label="butterworth filter")
plt.plot(w/ww, mavg_abs, label="moving average filter")


plt.title('Elliptical lowpass filter of order %d, decimation %d, ell_factor %d'%(N, decimation, ell_factor))
plt.xlabel('Frequency (MHz)')
plt.ylabel('Amplitude (dB)')
plt.grid(which='both', axis='both')
plt.fill([ws/ww*np.pi/decimation*ell_factor, max(w/ww), max(w/ww), ws*np.pi/ww/decimation*ell_factor], [max(h_abs), max(h_abs), -gstop, -gstop], '0.9', lw=0) # stop
plt.fill([wp/ww*np.pi/decimation*ell_factor, min(w/ww), min(w/ww), wp*np.pi/ww/decimation*ell_factor], [min(h_abs), min(h_abs), -gpass, -gpass], '0.9', lw=0) # stop
plt.axis([min(w/ww), max(w/ww), min(h_abs)-5, max(h_abs)+5])
plt.legend()
plt.show()
plt.savefig('c://lneuhaus//github//pyrpl//doc//specan_filter.png',DPI=300)

print "Final biquad coefficients [b0, b1, b2, a0, a1, a2]:"
for biquad in signal.zpk2sos(z, p, k):
    print biquad

Resulting filter

We see that a filter of 8th order, consisting of 4 sequential biquads is required. Since we do not require the span / sampling rate of the spectrum analyzer to be above roughly 5 MHz, we may implement the four biquads sequentially. Furthermore, for even lower values of the span, the filter can be fed with a reduced clock rate equal to the scope decimation factor divided by the variable ‘decimation’ in the filter design code above (4 in the example). For the aliasing of the lowpass filter passband not to cause problems in this case, we must in addition use the 4th order butterworth lowpass already available from the IQ module and the moving average filter of the scope. Then, as the plot shows, we can be sure that no aliasing occurs, given that no aliasing from the ADCs is present (should be guaranteed by analog Red Pitaya design).

The problem with our scheme is the complexity of introducing 2 (for the two quadratures) 4-fold biquads. This will not fit into the current design and must therefore be postponed to after the FPGA cleanup.

We could however opt for another temporary option, applicable only to stationary signals: Measure the spectrum twice or thrice with slightly shifted IQ demodulation frequency (at +- 10% of span and the actual center, as required above), and only plot the pointwise-minimum (with respect to the final frequency axis) of the obtained traces. This is simple and should be very effective (also to reduce the central peak at the demodulation freuqency), so i suggest we give it a try. Furthermore, it prepares the user that IQ spectra will only have 80% of the points in baseband mode, which will remain so after the implementation of the lowpass filter. The plot above shows that we do not have to worry about aliasing from multiple spans away if the bandwidth if the IQ module is se to the scope sampling rate (or slightly below). I am not aware that this method is used anywhere else, but do not see any serious problem with it.