Support (root) raised cosine filter impulse response.

applications/lofar2/model/pfb_os/dsp.py

@@ -117,6 +117,72 @@ def impulse_at_zero_crossing(x):
     return np.concatenate((np.array([0]), diff))
 
 
+###############################################################################
+# Window design
+###############################################################################
+
+def raised_cosine_response(Ntaps, Tsymbol, beta):
+    """Generate a raised cosine (RC) FIR filter impulse response.
+
+    Input:
+    . Ntaps : FIR filter length
+    . Tsymbol: symbol period in number of samples per symbol
+    . beta : Roll off factor in [0, 1.0], BW = (1 + beta) / Tsymbol, so:
+       - beta = 0.0: rectangular spectrum with BW = 1 / Tsymbol
+       - beta = 1.0: cosine spectrum with BW = 2 / Tsymbol
+    Return:
+    . hRc : impulse response of the raised cosine filter.
+    """
+    # time axis
+    tIndices = np.arange(Ntaps)
+    tCenter = (Ntaps - 1) / 2
+    t = tIndices - tCenter
+
+    # sinc term, can use array assignment because sinc(1 / 0) = 1
+    hRc = 1 / Tsymbol * np.sinc(t / Tsymbol)
+
+    # apply cos term, use for loop instead of array assignment, to detect divide by 0
+    for tI in tIndices:
+        t = tI - tCenter
+        if np.abs(t) != Tsymbol / (2 * beta):
+            hRc[tI] *= np.cos(np.pi * beta * t / Tsymbol) / (1 - (2 * beta * t / Tsymbol)**2)
+    return hRc
+
+
+def square_root_raised_cosine_response(Ntaps, Tsymbol, beta):
+    """Generate a square root raised cosine (RC) FIR filter impulse response.
+
+    Reference: [HARRIS section 4.3]
+    Input:
+    . Ntaps : FIR filter length
+    . Tsymbol: symbol period in number of samples per symbol
+    . beta : Roll off factor in [0, 1.0]
+    Return:
+    . hSrrc : impulse response of the square root raised cosine filter.
+    """
+    # time axis
+    tIndices = np.arange(Ntaps)
+    tCenter = (Ntaps - 1) / 2
+    t = tIndices - tCenter
+
+    # numerator term, using array assignment
+    hSrrc = 1 / Tsymbol * (4 * beta * t / Tsymbol * np.cos(np.pi * (1 + beta) * t / Tsymbol) +
+                                                    np.sin(np.pi * (1 - beta) * t / Tsymbol))
+
+    # apply denumerator term, use for loop instead of array assignment, to detect divide by 0
+    for tI in tIndices:
+        t = tI - tCenter
+        if t == 0.0:
+            hSrrc[tI] = 1 / Tsymbol * (1 + beta * (4 / np.pi - 1))
+        elif np.abs(t) == Tsymbol / (4 * beta):
+            hSrrc[tI] = 1 / Tsymbol * beta / np.sqrt(2) * \
+                ((1 + 2 / np.pi) * np.sin(np.pi / (4 * beta)) + \
+                 (1 - 2 / np.pi) * np.cos(np.pi / (4 * beta)))
+        else:
+            hSrrc[tI] /= (1 - (4 * beta * t / Tsymbol)**2) * (np.pi * t / Tsymbol)
+    return hSrrc
+
+
 ###############################################################################
 # FIR Filter design
 ###############################################################################
@@ -914,18 +980,23 @@ def resample(x, Nup, Ndown, coefs, verify=False):  # interpolate and decimate by
 # Plotting
 ###############################################################################
 
-def plot_time_response(h, name='', markers=False):
+def plot_time_response(h, title='', color='r', markers=False):
     """Plot time response (= impulse response, window, FIR filter coefficients).
 
     Input:
     . h: time response
+    . title: plot title
+    . color: curve color format character
     . markers: when True plot time sample markers in curve
     """
     if markers:
-        plt.plot(h, '-', h, 'o')
+        plt.plot(h, color + '-', h, color + 'o')
     else:
-        plt.plot(h, '-')
-    plt.title('Time response %s' % name)
+        plt.plot(h, color + '-')
+    if title:
+        plt.title(title)
+    else:
+        plt.title('Time response')
     plt.ylabel('Voltage')
     plt.xlabel('Sample')
     plt.grid(True)
@@ -1091,6 +1162,28 @@ def plot_power_spectrum(f, HF, fmt='r', fs=1.0, fLim=None, dbLim=None):
     plt.grid(True)
 
 
+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)
+    . 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
+    """
+    flabel = 'Frequency [fs = %f]' % fs
+
+    plt.plot(f * fs, np.abs(HF), fmt)
+    plt.title('Magnitude spectrum')
+    plt.ylabel('Voltage')
+    plt.xlabel(flabel)
+    if fLim:
+        plt.xlim(fLim)
+    if voltLim:
+        plt.ylim(voltLim)
+    plt.grid(True)
+
+
 def plot_two_power_spectra(f1, HF1, name1, f2, HF2, name2, fs=1.0, fLim=None, dbLim=None, showRoll=False):
     """Plot two power spectra in same plot for comparison