Easy way to implement a Root Raised Cosine (RRC) filter using Python & Numpy

Question

SciPy/Numpy seems to support many filters, but not the root-raised cosine filter. Is there a trick to easily create one rather than calculating the transfer function? An approximation would be fine as well.

Answer 1

The commpy package has several filters included with it. The order of return variables was switched in an earlier version (as of this edit, current version is 0.7.0). To install, foemphasized textllow instructions here or here.

Here's a use example for 1024 symbols of QAM16:

import numpy as np
from commpy.modulation import QAMModem
from commpy.filters import rrcosfilter
N = 1024  # Number of symbols
os = 8 #over sampling factor
# Create modulation. QAM16 makes 4 bits/symbol
mod1 = QAMModem(16)
# Generate the bit stream for N symbols
sB = np.random.randint(0, 2, N*mod1.num_bits_symbol)
# Generate N complex-integer valued symbols
sQ = mod1.modulate(sB)
sQ_upsampled = np.zeros(os*(len(sQ)-1)+1,dtype = np.complex64) 
sQ_upsampled[::os] = sQ
# Create a filter with limited bandwidth. Parameters:
#      N: Filter length in samples
#    0.8: Roll off factor alpha
#      1: Symbol period in time-units
#     24: Sample rate in 1/time-units
sPSF = rrcosfilter(N, alpha=0.8, Ts=1, Fs=over_sample)[1]
# Analog signal has N/2 leading and trailing near-zero samples
qW = np.convolve(sPSF, sQ_upsampled)

Here's some explanation of the parameters. N is the number of baud samples. You need 4 times as many bits (in the case of QAM) as samples. I made the sPSF array return with N elements so we can see the signal with leading and trailing samples. See the Wikipedia Root-raised-cosine filter page for explanation of parameter alpha. Ts is the symbol period in seconds and Fs is the number of filter samples per Ts. I like to pretend Ts=1 to keep things simple (unit symbol rate). Then Fs is the number of complex waveform samples per baud point.

If you use return element 0 from rrcosfilter to get the sample time indexes, you need to insert the correct symbol period and filter sample rate in Ts and Fs for the index values to be correctly scaled.

Answer 2

It would be nice to have the root-raised cosine filter standardized in a common package. Here is my implementation in the meantime based on commpy. It vectorized with numpy, and normalized without consideration of the symbol rate.

def raised_root_cosine(upsample, num_positive_lobes, alpha):
    """
    Root raised cosine (RRC) filter (FIR) impulse response.

    upsample: number of samples per symbol

    num_positive_lobes: number of positive overlaping symbols
    length of filter is 2 * num_positive_lobes + 1 samples

    alpha: roll-off factor
    """

    N = upsample * (num_positive_lobes * 2 + 1)
    t = (np.arange(N) - N / 2) / upsample

    # result vector
    h_rrc = np.zeros(t.size, dtype=np.float)

    # index for special cases
    sample_i = np.zeros(t.size, dtype=np.bool)

    # deal with special cases
    subi = t == 0
    sample_i = np.bitwise_or(sample_i, subi)
    h_rrc[subi] = 1.0 - alpha + (4 * alpha / np.pi)

    subi = np.abs(t) == 1 / (4 * alpha)
    sample_i = np.bitwise_or(sample_i, subi)
    h_rrc[subi] = (alpha / np.sqrt(2)) \
                * (((1 + 2 / np.pi) * (np.sin(np.pi / (4 * alpha))))
                + ((1 - 2 / np.pi) * (np.cos(np.pi / (4 * alpha)))))

    # base case
    sample_i = np.bitwise_not(sample_i)
    ti = t[sample_i]
    h_rrc[sample_i] = np.sin(np.pi * ti * (1 - alpha)) \
                    + 4 * alpha * ti * np.cos(np.pi * ti * (1 + alpha))
    h_rrc[sample_i] /= (np.pi * ti * (1 - (4 * alpha * ti) ** 2))

    return h_rrc

Answer 3

commpy doesn't seem to be released yet. But here is my nugget of knowledge.

beta = 0.20 # roll off factor

Tsample = 1.0 # sampling period, should at least twice the rate of the symbol

oversampling_rate = 8 # oversampling of the bit stream, this gives samples per symbol
# must be at least 2X the bit rate

Tsymbol = oversampling_rate * Tsample # pulse duration should be at least 2 * Ts
span = 50 # number of symbols to span, must be even
n = span*oversampling_rate # length of the filter = samples per symbol * symbol span

# t_step must be from -span/2 to +span/2 symbols.
# each symbol has 'sps' number of samples per second.
t_step = Tsample * np.linspace(-n/2,n/2,n+1) # n+1 to include 0 time

BW = (1 + beta) / Tsymbol
a = np.zeros_like(t_step)

for item in list(enumerate(t_step)):
    i,t = item 
    # t is n*Ts
    if (1-(2.0*beta*t/Tsymbol)**2) == 0:
        a[i] = np.pi/4 * np.sinc(t/Tsymbol)
        print 'i = %d' % i
    elif t == 0:
        a[i] = np.cos(beta * np.pi * t / Tsymbol)/ (1-(2.0*beta*t/Tsymbol)**2)
        print 't = 0 captured'
        print 'i = %d' % i 

    else:
        numerator = np.sinc( np.pi * t/Tsymbol )*np.cos( np.pi*beta*t/Tsymbol )
        denominator = (1.0 - (2.0*beta*t/Tsymbol)**2)
        a[i] =  numerator / denominator

#a = a/sum(a) # normalize total power

plot_filter = 0
if plot_filter == 1:

    w,h = signal.freqz(a)
    fig = plt.figure()
    plt.subplot(2,1,1)
    plt.title('Digital filter (raised cosine) frequency response')
    ax1 = fig.add_subplot(211)
    plt.plot(w/np.pi, 20*np.log10(abs(h)),'b')
    #plt.plot(w/np.pi, abs(h),'b')
    plt.ylabel('Amplitude (dB)', color = 'b')
    plt.xlabel(r'Normalized Frequency ($\pi$ rad/sample)')

    ax2 = ax1.twinx()
    angles = np.unwrap(np.angle(h))
    plt.plot(w/np.pi, angles, 'g')
    plt.ylabel('Angle (radians)', color = 'g')
    plt.grid()
    plt.axis('tight')
    plt.show()


    plt.subplot(2,1,2)
    plt.stem(a)
    plt.show()