SciPy/Numpy 似乎支持许多过滤器,但不支持根升余弦过滤器。有没有一种技巧可以轻松创建一个而不是计算传递函数?近似值也可以。
commpy
包中包含多个过滤器。返回变量的顺序在早期版本中已更改(截至本次编辑,当前版本为 0.8.0)。要安装,请按照说明此处或此处。
这是 QAM16 的 1024 个符号的示例:
import numpy as np
from commpy.modulation import QAMModem
from commpy.filters import rrcosfilter
N = 1024 # Number of symbols. Also, the filter length in samples.
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
# os: Sample rate in 1/time-units
sPSF = rrcosfilter(N, alpha=0.8, Ts=1, Fs=os)[1]
# Analog signal has N/2 leading and trailing near-zero samples
qW = np.convolve(sPSF, sQ_upsampled)
这里是一些参数的解释。
N
是波特率样本数。您需要 4 倍于样本的位数(在 QAM 的情况下)。我用 sPSF
元素返回 N
数组,这样我们就可以看到带有前导和尾随样本的信号。有关参数 alpha
的说明,请参阅维基百科根升余弦滤波器页面。
Ts
是以秒为单位的符号周期,Fs
是每个 Ts
的滤波器样本数。我喜欢假装 Ts=1
让事情变得简单(单位符号率)。那么 Fs
是每个波特点的复杂波形样本数。
如果使用
rrcosfilter
中的返回元素 0 来获取采样时间索引,则需要在 Ts
和 Fs
中插入正确的符号周期和滤波器采样率,以便正确缩放索引值。
如果能在通用包中标准化根升余弦滤波器,那就太好了。这是我同时基于 comppy 的实现。它使用 numpy 进行矢量化,并在不考虑符号率的情况下进行归一化。
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
commpy 好像还没有发布。但这是我的知识点。
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()
我认为正确的反应是产生期望的脉冲响应。对于升余弦滤波器,函数为
h(n) = (sinc(n/T)*cos(pi * alpha* n /T)) / (1-4*(alpha*n/T)**2)
选择过滤器的点数并生成权重。
output = scipy.signal.convolve(signal_in, h)
这与 CommPy 中的函数基本相同,但代码要小得多:
def rcosfilter(N, beta, Ts, Fs):
t = (np.arange(N) - N / 2) / Fs
return np.where(np.abs(2*t) == Ts / beta,
np.pi / 4 * np.sinc(t/Ts),
np.sinc(t/Ts) * np.cos(np.pi*beta*t/Ts) / (1 - (2*beta*t/Ts) ** 2))
SciPy 将支持任何过滤器。只需计算脉冲响应并使用任何适当的 scipy.signal 滤波器/卷积函数即可。