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.

Solution 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.

Solution 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

Solution 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()

Comments

[Dan Sandberg]

Great, thanks! A description of the parameters is available here: commpy.readthedocs.org/en/latest/generated/…

[nedlrichards]

I don't see how the rrc filter can be implemented with a sinc function, based on wikipedia definition.

[Yarden Cohen]

can you please explain those lengthes and parameters? it looks like you are the only usage example for that library around.. what is mod1.num_bits_symbol*N*M/4? and why rrcosfilter() gets N*4? also how does the Ts=1 and Fs=24 affect here?

[Sergey]

@nedlrichards found this equation here ntrs.nasa.gov/api/citations/20120008631/downloads/…

[nedlrichards]

Thanks for the reference! I think that this formulation is a form of Eq. 2 in the pdf, which is a raised cosine filter. This clarified for me what the original answer was doing. However, since the original question concerned the raised root cosine filter, I do not think this is the requested solution.

[Okan Erturk]

This is Raised Cosine, not root-Raised Cosine!