7

SciPy / Numpyは多くのフィルターをサポートしているようですが、ルートレイズドコサインフィルターはサポートしていません。伝達関数を計算するのではなく、簡単に作成するためのトリックはありますか？概算も問題ありません。

4

6 に答える 6

7

commpyパッケージにはいくつかのフィルターが含まれています。以前のバージョンでは、戻り変数の順序が入れ替わっていました（この編集の時点では、現在のバージョンは0.7.0です）。インストールするには、ここまたはここでテキストの説明を強調しました。

QAM16の1024シンボルの使用例を次に示します。

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)

パラメータの説明は次のとおりです。Nボーサンプルの数です。サンプルの4倍のビット（QAMの場合）が必要です。sPSF配列を要素で返すようにしたNので、先頭と末尾のサンプルで信号を確認できます。パラメーターの説明については、WikipediaのRoot-raised-cosineフィルターページalphaを参照してください。Tsは秒単位のシンボル周期であり、Fsはあたりのフィルターサンプルの数Tsです。Ts=1私は物事を単純に保つふりをするのが好きです（単位シンボルレート）。次にFs、ボーポイントあたりの複素波形サンプルの数です。

からの戻り要素を使用rrcosfilterしてサンプル時間インデックスを取得する場合は、正しいシンボル期間を挿入し、インデックス値を正しくスケーリングするためにサンプルレートTsをフィルタリングする必要があります。Fs

于 2014-09-15T22:44:26.280 に答える
2

ルートレイズドコサインフィルターを共通のパッケージで標準化しておくと便利です。これが、commpyに基づいた当面の私の実装です。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
于 2020-05-26T21:09:07.250 に答える
1

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()
于 2015-03-09T20:34:20.513 に答える
0

正しい応答は、欲求インパルス応答を生成することだと思います。レイズドコサインフィルターの場合、関数は次のようになります。

h(n) = (sinc(n/T)*cos(pi * alpha* n /T)) / (1-4*(alpha*n/T)**2)

フィルタのポイント数を選択し、重みを生成します。

output = scipy.signal.convolve(signal_in, h)
于 2015-03-03T22:21:34.897 に答える
0

これは基本的に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))
于 2019-06-10T17:13:38.150 に答える
-2

SciPyはすべてのフィルターをサポートします。インパルス応答を計算し、適切なscipy.signalフィルター/畳み込み関数のいずれかを使用するだけです。

于 2013-04-03T09:26:55.337 に答える