CommPy
CommPy copied to clipboard
rrcosfilter with Fs=3
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()

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