Improve frequency axis.

applications/lofar2/model/rtdsp/fourier.py

@@ -110,7 +110,7 @@ def dtft(coefs, Ndtft=None, zeroCenter=True, fftShift=True):
       else use 0 to 1.0.
     Return:
     . h: zero padded coefs
-    . f: normalized frequency axis for HF, fs = 1
+    . fn: normalized frequency axis for HF, fs = 1
     . HF: dtft(h), the frequency transfer function of h
     """
     c_interpolate = 10
@@ -125,12 +125,12 @@ def dtft(coefs, Ndtft=None, zeroCenter=True, fftShift=True):
     # DFT
     HF = np.fft.fft(h)
     # Normalized frequency axis, fs = 1, ws = 2 pi
-    f = np.arange(0, 1, 1 / Ndtft)  # f = 0,pos, neg
+    fn = np.arange(0, 1, 1 / Ndtft)  # fn = 0,pos, neg
     if fftShift:
-        # FFT shift to center HF, f = neg, 0,pos
-        f = f - 0.5
+        # FFT shift to center HF, fn = neg, 0,pos
+        fn = fn - 0.5
         HF = np.roll(HF, Ndtft // 2)
-    return h, f, HF
+    return h, fn, HF
 
 
 ###############################################################################
@@ -141,7 +141,8 @@ def estimate_gain_at_frequency(f, HF, freq):
     """Find gain = abs(HF) at frequency in f that is closest to freq.
 
     Input:
-    . f, HF: frequency axis and spectrum as from dtft() with fftShift=True
+    . f: frequency axis as f = fn * fs for fn from dtft() with fftShift=True
+    . HF: spectrum as from dtft() with fftShift=True
     . freq : frequency to look for in f
     Return:
     . fIndex: index of frequency in f that is closest to freq
@@ -158,19 +159,26 @@ def estimate_gain_at_frequency(f, HF, freq):
 def estimate_frequency_at_gain(f, HF, gain):
     """Find positive frequency in f that has gain = abs(HF) closest to gain.
 
+    Look only in positive frequencies and from highest frequency to lowest
+    frequency to avoid band pass ripple gain variation in LPF.
+
     Input:
-    . f, HF: frequency axis and spectrum as from dtft() with fftShift=True
+    . f: frequency axis as f = fn * fs for fn from dtft() with fftShift=True
+    . HF: spectrum as from dtft() with fftShift=True
     . gain : gain to look for in HF
     Return:
     . fIndex: index of frequency in f that has abs(HF) closest to gain
     . fValue: frequency at fIndex
     . fGain : gain = abs(HF) at fIndex
     """
+    # Look only in positive frequencies
     pos = len(f[f <= 0])
     HFpos = HF[pos:]
     HFgain = np.abs(HFpos)
     HFdiff = np.abs(HFgain - gain)
-    fIndex = pos + HFdiff.argmin()
+    # Flip to avoid detections in LPF band pass
+    HFdiff = np.flipud(HFdiff)
+    fIndex = pos + len(HFpos) - 1 - HFdiff.argmin()
     fValue = f[fIndex]
     fGain = np.abs(HF[fIndex])
     return fIndex, fValue, fGain

applications/lofar2/model/rtdsp/plotting.py

@@ -57,7 +57,7 @@ def plot_time_response(h, title='', color='r', markers=False):
     plt.grid(True)
 
 
-def plot_iir_filter_analysis(b, a, fs=1, whole=False, Ntime=100, step=False, log=False, show=[]):
+def plot_iir_filter_analysis(b, a, fs=1.0, whole=False, Ntime=100, step=False, log=False, show=[]):
     """Plot and print iir filter analysis results.
 
     Input:
@@ -137,25 +137,22 @@ def plot_iir_filter_analysis(b, a, fs=1, whole=False, Ntime=100, step=False, log
     return z, p, k
 
 
-def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
+def plot_spectra(fn, HF, fs=1.0, fLim=None, dbLim=None):
     """Plot spectra for power, magnitude, phase, real, imag
 
     Input:
-    . f: normalized frequency axis for HF (fs = 1)
+    . fn: normalized frequency axis for HF (fs = 1)
     . HF: spectrum, e.g. frequency transfer function HF = DTFT(h)
-    . fs: sample frequency in Hz, scale f by fs, fs >= 1
+    . fs: sample frequency in Hz
     """
     Hmag = np.abs(HF)
     Hphs = np.unwrap(np.angle(HF))
     Hpow_dB = pow_db(HF)  # power response
-    fn = f * fs
-    if fs > 1:
-        flabel = 'Frequency [fs / %d]' % fs
-    else:
-        flabel = 'Frequency [fs]'
+    f = fn * fs  # scale fn by fs
+    flabel = 'Frequency [fs = %f]' % fs
 
     plt.figure(1)
-    plt.plot(fn, Hpow_dB)
+    plt.plot(f, Hpow_dB)
     plt.title('Power spectrum')
     plt.ylabel('Power [dB]')
     plt.xlabel(flabel)
@@ -166,8 +163,8 @@ def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
     plt.grid(True)
 
     plt.figure(2)
-    plt.plot(fn, HF.real, 'r')
-    plt.plot(fn, HF.imag, 'g')
+    plt.plot(f, HF.real, 'r')
+    plt.plot(f, HF.imag, 'g')
     plt.title('Complex voltage spectrum')
     plt.ylabel('Voltage')
     plt.xlabel(flabel)
@@ -177,7 +174,7 @@ def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
     plt.grid(True)
 
     plt.figure(3)
-    plt.plot(fn, Hmag)
+    plt.plot(f, Hmag)
     plt.title('Magnitude spectrum')  # = amplitude
     plt.ylabel('Voltage')
     plt.xlabel(flabel)
@@ -186,7 +183,7 @@ def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
     plt.grid(True)
 
     plt.figure(4)
-    plt.plot(fn, Hphs)
+    plt.plot(f, Hphs)
     plt.title('Phase spectrum (note -1: pi = -pi)')
     plt.ylabel('Phase [rad]')
     plt.xlabel(flabel)
@@ -195,18 +192,19 @@ def plot_spectra(f, HF, fs=1.0, fLim=None, dbLim=None):
     plt.grid(True)
 
 
-def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None):
+def plot_power_spectrum(fn, HF, fmt='r', fs=1.0, fLim=None, dbLim=None):
     """Plot power spectrum
 
     Input:
-    . f: normalized frequency axis for HF (fs = 1)
+    . fn: normalized frequency axis for HF (fs = 1)
     . HF: spectrum, e.g. frequency transfer function HF = DTFT(h)
     . fmt: curve format string
-    . fs: sample frequency in Hz, scale f by fs, fs >= 1
+    . fs: sample frequency in Hz
     """
+    f = fn * fs  # scale fn by fs
     flabel = 'Frequency [fs = %f]' % fs
 
-    plt.plot(f * fs, pow_db(HF), fmt)
+    plt.plot(f, pow_db(HF), fmt)
     plt.title('Power spectrum')
     plt.ylabel('Power [dB]')
     plt.xlabel(flabel)
@@ -221,11 +219,12 @@ def plot_magnitude_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, voltLim=None):
     """Plot magnitude (= voltage) spectrum
 
     Input:
-    . f: normalized frequency axis for HF (fs = 1)
+    . fn: normalized frequency axis for HF (fs = 1)
     . HF: spectrum, e.g. frequency transfer function HF = DTFT(h)
     . fmt: curve format string
-    . fs: sample frequency in Hz, scale f by fs, fs >= 1
+    . fs: sample frequency in Hz
     """
+    f = fn * fs  # scale fn by fs
     flabel = 'Frequency [fs = %f]' % fs
 
     plt.plot(f * fs, np.abs(HF), fmt)
@@ -239,28 +238,27 @@ def plot_magnitude_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, voltLim=None):
     plt.grid(True)
 
 
-def plot_two_power_spectra(f1, HF1, name1, f2, HF2, name2, fs=1.0, fLim=None, dbLim=None, showRoll=False):
+def plot_two_power_spectra(fn1, HF1, name1, fn2, HF2, name2, fs=1.0, fLim=None, dbLim=None, showRoll=False):
     """Plot two power spectra in same plot for comparison
 
     Input:
-    . f1,f2: normalized frequency axis for HF1, HF2 (fs = 1)
+    . fn1,fn2: normalized frequency axis for HF1, HF2 (fs = 1)
     . HF1, HF2: spectrum, e.g. frequency transfer function HF = DTFT(h)
-    . fs: sample frequency in Hz, scale f by fs, fs >= 1
+    . fs: sample frequency in Hz
     """
-    if fs > 1:
-        flabel = 'Frequency [fs / %d]' % fs
-    else:
-        flabel = 'Frequency [fs]'
+    f1 = fn1 * fs  # scale fn1 by fs
+    f2 = fn2 * fs  # scale fn1 by fs
+    flabel = 'Frequency [fs = %f]' % fs
 
     if showRoll:
-        plt.plot(f1 * fs, pow_db(HF1), 'r',
-                 f1 * fs, np.roll(pow_db(HF1), len(f1) // 2), 'r--',
-                 f2 * fs, pow_db(HF2), 'b',
-                 f2 * fs, np.roll(pow_db(HF2), len(f2) // 2), 'b--')
+        plt.plot(f1, pow_db(HF1), 'r',
+                 f1, np.roll(pow_db(HF1), len(f1) // 2), 'r--',
+                 f2, pow_db(HF2), 'b',
+                 f2, np.roll(pow_db(HF2), len(f2) // 2), 'b--')
         plt.legend([name1, '', name2, ''])
     else:
-        plt.plot(f1 * fs, pow_db(HF1), 'r',
-                 f2 * fs, pow_db(HF2), 'b')
+        plt.plot(f1, pow_db(HF1), 'r',
+                 f2, pow_db(HF2), 'b')
         plt.legend([name1, name2])
     plt.title('Power spectrum')
     plt.ylabel('Power [dB]')