CommPy icon indicating copy to clipboard operation
CommPy copied to clipboard

rrcosfilter with Fs=3

Open philn16 opened this issue 6 years ago • 1 comments
trafficstars

I've noticed a bizzar frequency response from the rrcosfilter. I suspect this is an edge case due to the if/else in the module. The time waveform doesn't look unusual.

import matplotlib.pyplot as plt
import scipy.fftpack as fftpack
from commpy import rrcosfilter
import numpy as np
import scipy

def padded_centered_fft(taps,mintaps=1024,db=True,db_min=-300,fs=1):
	""" returns the fft of the zero padded taps, and centers it """
	if len(taps) < mintaps:
		taps=np.concatenate((taps,np.zeros(mintaps-len(taps))))
	dfreq, fft = fftpack.fftshift(fftpack.fftfreq(len(taps))),fftpack.fftshift(scipy.fft(taps))
	fft_abs=np.abs(fft)
	fft_db = 20*np.log10(fft_abs+max(np.abs(fft))*(10**((db_min-50)/20)) )
	fft_db[np.where(fft_db < db_min)]=db_min
	dfreq*=fs
	return (dfreq, fft_db) if db else (dfreq,fft)

plt.plot(*padded_centered_fft(rrcosfilter(N=400,alpha=0.15,Ts=1,Fs=3.01)[1]),label='3.01')
plt.plot(*padded_centered_fft(rrcosfilter(N=400,alpha=0.15,Ts=1,Fs=3)[1]),label='3')
plt.plot(*padded_centered_fft(rrcosfilter(N=400,alpha=0.15,Ts=1,Fs=2.99)[1]),':',label='2.99')
plt.legend(loc='best')
plt.show()

response

philn16 avatar Aug 02 '19 19:08 philn16

here's the rrcosfilter function from the commpy I'm using:

def rrcosfilter(N, alpha, Ts, Fs):
    """
    Generates a root raised cosine (RRC) filter (FIR) impulse response.

    Parameters
    ----------
    N : int
        Length of the filter in samples.

    alpha : float
        Roll off factor (Valid values are [0, 1]).

    Ts : float
        Symbol period in seconds.

    Fs : float
        Sampling Rate in Hz.

    Returns
    ---------

    time_idx : 1-D ndarray of floats
        Array containing the time indices, in seconds, for
        the impulse response.

    h_rrc : 1-D ndarray of floats
        Impulse response of the root raised cosine filter.
    """

    T_delta = 1/float(Fs)
    time_idx = ((np.arange(N)-N/2))*T_delta
    sample_num = np.arange(N)
    h_rrc = np.zeros(N, dtype=float)

    for x in sample_num:
        t = (x-N/2)*T_delta
        if t == 0.0:
            h_rrc[x] = 1.0 - alpha + (4*alpha/np.pi)
        elif alpha != 0 and t == Ts/(4*alpha):
            h_rrc[x] = (alpha/np.sqrt(2))*(((1+2/np.pi)* \
                    (np.sin(np.pi/(4*alpha)))) + ((1-2/np.pi)*(np.cos(np.pi/(4*alpha)))))
        elif alpha != 0 and t == -Ts/(4*alpha):
            h_rrc[x] = (alpha/np.sqrt(2))*(((1+2/np.pi)* \
                    (np.sin(np.pi/(4*alpha)))) + ((1-2/np.pi)*(np.cos(np.pi/(4*alpha)))))
        else:
            h_rrc[x] = (np.sin(np.pi*t*(1-alpha)/Ts) +  \
                    4*alpha*(t/Ts)*np.cos(np.pi*t*(1+alpha)/Ts))/ \
                    (np.pi*t*(1-(4*alpha*t/Ts)*(4*alpha*t/Ts))/Ts)

    return time_idx, h_rrc

philn16 avatar Aug 02 '19 19:08 philn16