Spectral inversion.

libraries/dsp/si/src/python/si.py

+#! /usr/bin/env python3
+###############################################################################
+#
+# Copyright 2022
+# ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+###############################################################################
+
+# Author: Eric Kooistra
+# Date: Nov 2022
+# Purpose:
+#   Derive RF frequency of subband center dependend on subband index, Nyquist
+#   zone and spectral inversion.
+# Description:
+#   The > 1 st Nyquist zones are digitized using sub sampling. The analog
+#   filter in the receiver blocks the other Nyquist zones, so that only the
+#   wanted Nyquist zone is digitized.
+#   The sub sampling act as a mixer that shifts all Nyquist zones down to the
+#   first Nyquist zone. The sub sampling causes that the frequency band for
+#   all even Nyquist zones gets flipped from hi to lo. This frequency flip can
+#   be undone by enabling the spectral inversion for those even Nyquist zones.
+#   For ADC sample frequenc f_sample = 200 MHz:
+#     1st Nyquist zone: f =   0 - 100 MHz
+#     2nd Nyquist zone: f = 100 - 200 MHz
+#     3rd Nyquist zone: f = 200 - 300 MHz
+
+import argparse
+import textwrap
+
+import numpy as np
+import matplotlib
+matplotlib.use('tkagg')
+import matplotlib.pyplot as plt
+
+# Parse arguments to derive user parameters
+_parser = argparse.ArgumentParser(
+    description="".join(textwrap.dedent("""\
+        Calculate RF center frequency for subbands dependend on the Nyquist
+        zone and the spectral inversion control.
+
+        Apply spectral inversion for 2nd Nyquist zone to have increasing
+        subband index correspond to increasing RF frequency:
+
+        > python si.py --si 0 --zi 0 -N 16
+        > python si.py --si 1 --zi 1 -N 16
+        > python si.py --si 0 --zi 2 -N 16
+
+        \n""")),
+    formatter_class=argparse.RawTextHelpFormatter)
+_parser.add_argument('--si', default=0, type=int, help='Spectral inversion control 0 = keep band, 1 = flip band')
+_parser.add_argument('--zi', default=0, type=int, help='Nyquist zone index 0 = 1st, 1 = 2nd, 2 = 3rd, etc')
+_parser.add_argument('-N', default=512, type=int, help='N_sub')
+args = _parser.parse_args()
+
+spectral_inv = args.si
+nyquist_zone_index = args.zi
+N_sub = args.N
+
+subbands = np.arange(N_sub)   # subband indices
+
+# SDP parameters
+N_complex = 2                 # Two complex parts, real and imag
+N_fft = N_sub * N_complex     # Number of points of subband filterbank FFT
+f_sample = 200e6              # Hz (ADC sample frequency)
+BW_RF = f_sample / N_complex  # = 100 MHz, sampled RF band width of each
+                              # Nyquist zone,
+f_sub = f_sample / N_fft      # = 195312.5 Hz, subband frequency
+
+# Determine netto result of sub sampling and spectral inversion control
+zone_inv = nyquist_zone_index % N_complex
+
+def boolean_xor(a, b):
+    return (a and not b) or (b and not a)
+
+def natural_xor(a, b):
+    if boolean_xor(a, b):
+        return 1
+    else:
+        return 0
+
+bw_inv = natural_xor(spectral_inv, zone_inv)
+
+# Subband center RF frequency
+f_lo = nyquist_zone_index * BW_RF
+
+n = subbands
+if bw_inv:
+    n = N_sub - 1 - subbands
+
+f_sub_rf = f_lo + n * f_sub
+
+# Plot results
+figNr = 0
+
+figNr += 1
+plt.figure(figNr)
+plt.plot(subbands, f_sub_rf / 1e6, 'o')
+plt.title("spectral_inv = %d, nyquist_zone_index = %d" % (spectral_inv, nyquist_zone_index))
+plt.xlabel("Subband index (range 0:%d)" % (N_sub - 1))
+plt.ylabel("Subband RF frequency [MHz]")
+plt.grid()
+
+plt.show()