【Python】commpyとFuncAnimationでQPSK変調のIQ信号波形を動画で描いてみる

今回はPythonのFuncAnimationを使って↓な動画を作成します。
f:id:taekwongineer:20210619231754g:plain

QPSK変調信号のレイズドコサインフィルタでの帯域制限前後のIch、Qchの信号の時間波形とIch、Qchの散布図です。

QPSK変調信号を作成するのに便利なcommpyとMatplotlibのグラフで動画を作成するFuncAnimationの使い方を備忘録としてまとめます。

●この記事の目次

1. まずはソースコード
2. commpyを活用した信号処理
- 2.1 commpy.utilities.upsample
- 2.2 commpy.filters.rcosfilter
3. FuncAnimationによる動画の作成
- 3.1 動画の作成方法
- 3.2 動画の保存方法
4. まとめ

1. まずはソースコード

いきなりですがソースコードはこちら

import numpy as np
import commpy as cp
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation as fa

######## ①QPSK変調信号生成 ########
####　変数定義 ####
N = 10000                   # 送信ビット数
k = 8                       # アップサンプリングの倍率
SymbolRate = 1000           # シンボルレート[Hz]
SampleRate = k * SymbolRate # サンプリングレート[Hz]
Tap = k * 4                 # レイズドコサインフィルタのタップ数
Alpha = 0.5                 # ロールオフ率

#### 送信ビットを生成 ####
TxData = np.random.randint(0, 2, N)

#### QPSKのシンボルマッピング ####
# 00 : (+1+i)/√2
# 01 : (+1-i)/√2
# 10 : (-1+i)/√2
# 11 : (-i-i)/√2
QpskSym = np.array([])    # QPSK変調のシンボルを格納する配列
tmp = 0  # QPSK変調のシンボルの計算結果を格納する一時変数
for i in range(0, N, 2):
    
    if TxData[i] == 0:
        tmp = 1
    else:
        tmp = -1
        
    if TxData[i+1] == 0:
        tmp = tmp + 1j
    else:
        tmp = tmp - 1j
    
    tmp = tmp / np.sqrt(2)
    
    QpskSym = np.append(QpskSym, tmp)

#### アップサンプリング ####
UpSamp = cp.utilities.upsample(QpskSym, k)

#### レイズドコサインフィルタによる帯域制限　####
# レイズドコサインフィルタを作成
t, Filter = cp.filters.rcosfilter(Tap, Alpha, 1/SymbolRate, SampleRate)

# アップサンプリングされた信号とレイズドコサインフィルタを畳み込み積分
FilteredSamp = np.convolve(UpSamp, Filter)


######## ②QPSK変調のIQ信号を可視化 ########
#### 変数定義 ####
L = Tap * 2                 # グラフに表示する横軸のサンプル数
n = 0                       # 現在のサンプル

#　グラフ表示用の配列の作成
IchPlotIn = np.zeros(0)     # Ichのレイズドコサインフィルタ入力
QchPlotIn = np.zeros(0)     # Qchのレイズドコサインフィルタ入力
IchPlotOut = np.zeros(0)    # Ichのレイズドコサインフィルタ出力
QchPlotOut = np.zeros(0)    # Qchのレイズドコサインフィルタ出力
TimePlot = np.zeros(0)      # 時刻

Ts = np.arange(0,len(FilteredSamp)) / SampleRate * 1000 #　時刻[ms]

# グラフ描画エリアの作成
fig = plt.figure(figsize=[12, 8])
plt.subplots_adjust(wspace=0.3, hspace=0.8)

#### FuncAnimation関数に入力する表示更新処理用関数 ####
def update(f):
    global IchPlotIn, QchPlotIn, IchPlotOut, QchPlotOut, TimePlot, n, FilteredSamp
    
    # グラフの表示内容をクリア
    plt.clf()
    
    # グラフの描画エリアを定義
    ax_IchIn = fig.add_subplot(421)
    ax_IchIn.set_ylim(-1.5, 1.5)
    ax_IchIn.set_title("Ich Signal(Input)", fontsize=10)
    ax_IchIn.set_xlabel("Time[ms]", fontsize=8)
    ax_IchIn.set_ylabel("Level", fontsize=8)
    
    ax_IchOut = fig.add_subplot(423)
    ax_IchOut.set_ylim(-1.5, 1.5)
    ax_IchOut.set_title("Ich Signal(Output)", fontsize=10)
    ax_IchOut.set_xlabel("Time[ms]", fontsize=8)
    ax_IchOut.set_ylabel("Level", fontsize=8)
    
    ax_QchIn = fig.add_subplot(425)
    ax_QchIn.set_ylim(-1.5, 1.5)
    ax_QchIn.set_title("Qch Signal(Input)", fontsize=10)
    ax_QchIn.set_xlabel("Time[ms]", fontsize=8)
    ax_QchIn.set_ylabel("Level", fontsize=8)
    
    ax_QchOut = fig.add_subplot(427)
    ax_QchOut.set_ylim(-1.5, 1.5)
    ax_QchOut.set_title("Qch Signal(Output)", fontsize=10)
    ax_QchOut.set_xlabel("Time[ms]", fontsize=8)
    ax_QchOut.set_ylabel("Level", fontsize=8)
    
    ax_IQ = fig.add_subplot(122, aspect=1) 
    ax_IQ.set_xlim(-1.5, 1.5)
    ax_IQ.set_ylim(-1.5, 1.5)
    ax_IQ.set_title("IQ Trajectory", fontsize=10)
    ax_IQ.set_xlabel("Ich Level", fontsize=8)
    ax_IQ.set_ylabel("Qch Level", fontsize=8)
    
    # グラフ表示用の配列の要素数に応じて処理を分岐
    if len(IchPlotIn) < L:
        # 要素数がLに満たない場合は、X軸の表示範囲は固定
        ax_IchIn.set_xlim(0, Ts[L-1])
        ax_IchOut.set_xlim(0, Ts[L-1])
        ax_QchIn.set_xlim(0, Ts[L-1])
        ax_QchOut.set_xlim(0, Ts[L-1])
        
        # 要素数がLに到達するまで、配列に要素を追加
        IchPlotIn = np.append(IchPlotIn, UpSamp[n].real)
        IchPlotOut = np.append(IchPlotOut, FilteredSamp[n].real)
        QchPlotIn = np.append(QchPlotIn, UpSamp[n].imag)
        QchPlotOut = np.append(QchPlotOut, FilteredSamp[n].imag)
        
        TimePlot = np.append(TimePlot, Ts[n])
        
    else:
        # 配列の要素数がLに到達した後は、先頭の要素を削除後に、次の要素を追加
        IchPlotIn = IchPlotIn[1:L]
        IchPlotIn = np.append(IchPlotIn, UpSamp[n].real)
        
        IchPlotOut = IchPlotOut[1:L]
        IchPlotOut = np.append(IchPlotOut, FilteredSamp[n].real)
        
        QchPlotIn = QchPlotIn[1:L]
        QchPlotIn = np.append(QchPlotIn, UpSamp[n].imag)
        
        QchPlotOut = QchPlotOut[1:L]
        QchPlotOut = np.append(QchPlotOut, FilteredSamp[n].imag)
        
        TimePlot = TimePlot[1:L]
        TimePlot = np.append(TimePlot, Ts[n])
        
    ax_IchIn.plot(TimePlot, IchPlotIn, "-xb")
    ax_IchIn.grid(True)
    
    ax_IchOut.plot(TimePlot, IchPlotOut, "-xm")
    ax_IchOut.grid(True)
    
    ax_QchIn.plot(TimePlot, QchPlotIn, "-xb")
    ax_QchIn.grid(True)
    
    ax_QchOut.plot(TimePlot, QchPlotOut, "-xm")
    ax_QchOut.grid(True)
    
    ax_IQ.plot(IchPlotOut, QchPlotOut)
    ax_IQ.plot(np.array([1,-1,1,-1])/np.sqrt(2), np.array([1,1,-1,-1])/np.sqrt(2), "ro")
    ax_IQ.plot(FilteredSamp[n].real, FilteredSamp[n].imag, "xk")
    ax_IQ.grid(True)
    
    # 参照する配列のインデックスを更新
    n = n + 1

#### FuncAnimationによるアニメーション描画の実行 ####
ani = fa(fig, update, interval=10, frames=500)
#ani.save("output.gif", writer="imagemagick")
plt.show()

2. commpyを活用した信号処理

commpyはMATLABのCommunications Toolboxのような機能を持ったパッケージで、

・変復／復調
・誤り訂正符号化／誤り訂正
・フェージング生成
・フィルタ生成

などなど、ディジタル無線通信のシミュレーションで多用する各種機能の関数群が含まれています。

commpy.readthedocs.io

MATLABもCommunications Toolboxも個人のお小遣いで買えるような代物ではありませんからね、このようなパッケージがあるのは非常にありがたいことです。

今回はこれらの中から
・commpy.utilities.upsample　：アップサンプリング
・commpy.filters.rcosfilter　　：レイズドコサインフィルター
の機能を使用します。

2.1 commpy.utilities.upsample

upsampleは、その名の通りアップサンプリングを行う関数です。

使い方は

Output = commpy.utilities.upsample( Input, N)

・引数　　　Input (1-D ndarray(float)) ：入力サンプル
　　　　　　　 N (int)　：アップサンプリングの倍率
・戻り値　Output (1-D ndarray(float)) ：アップサンプリング後のデータ列

となります。

例えば、Input=[1, -1, 1]の配列を、N=4でアップサンプリングすると、

Output = [1, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0]

のようになります。

今回のサンプルのコードでは、42行目で
UpSamp = cp.utilities.upsample(QpskSym, k)
の形で使用しています。

2.2 commpy.filters.rcosfilter

rcosfilterはレイズドコサインフィルターを生成する関数です。

使い方は

time_idx, rc_filter = commpy.filters.rcosfilter(N, Alpha, Ts, Fs)

・引数　　　 N(int)　：生成するレイズドコサインフィルターのタップ数
　　　　Alpha(float) 　：ロールオフ率（0.0～1.0）
　　　　　　Ts(float)　： 1シンボル時間[s]（シンボルレート[Hz]の逆数）
　　　　　　Fs(float)　：サンプリングレート[Hz]
・戻り値　　time_idx(1-D ndarray(float)) : レイズドコサインフィルタの時刻インデックス
　　　　　　rc_filter(1-D ndarray(float)) : レイズドコサインフィルタのインパルス応答

となります。

例えば
・N = 48
・Alpha = 0.5
・Ts = 0.001(=1/1000)
・Fs = 8000
とすると、出力されるレイズドコサインフィルタのインパルス応答は
f:id:taekwongineer:20210620221041p:plain
のようになります。

グラフの横軸は、戻り値のtime_idxを使用しています。

今回のサンプルのコードでは、46行目で
t, Filter = cp.filters.rcosfilter(Tap, Alpha, 1/SymbolRate, SampleRate)
の形で使用しています。

3. FuncAnimationによる動画の作成

次に、2.で作成したQPSK変調信号を動画として描画する方法について説明します。

今回は、matplotlib.animation.FuncAnimatitonを使用し、動画を作成します。

3.1 動画の作成方法

ここでは、FuncAnimationの簡単な使用方法のみ示します。

サンプルコードでは、

ani = fa(fig, update, interval=10, frames=500)
・引数：　fig(figure)　：　動画の描画対象となるグラフエリアのfigureオブジェクト
　　　　　update(function)　：　グラフの描画の更新処理を行う関数
　　　　　interval　：　描画の更新周期(ms)（指定がない場合のデフォルト値は200ms）
　　　　　frames　：　動画のフレーム数（指定がない場合は無限）

の形式で使用しています。

他にもいろいろ引数で条件は指定できるのですが、詳細は↓を参照してください。
matplotlib.org

3.2 動画の保存方法

動画を保存する場合は、

FuncAnimation.save(Filename, writer)
・引数：　Filename（string)　：　保存ファイル名（gif)
　　　　　writer（string）　：　gifファイルを作成するwriter

と記述します。

今回のサンプルコードでは、wriiterにimagemagickを使用していますが、デフォルトではインストールされていないので、別途インストールおよび設定が必要となります。
imagemagickのインストールと設定に関しては、↓のサイトに分かりやすく説明が書かれています。
imagingsolution.net

4. まとめ

今回はcommpyとFuncAnimationを使用して、QPSK変調波形を動画で描く方法についてまとめました。
commpyについては、誤り訂正符号化／誤り訂正の使い方を抑えておくと、誤り訂正付きのBERのシミュレーションを作成するのがとても楽になりそうなので、使い方を調査していきたいと思います。

無線通信エンジニアの備忘録

無線通信だったり、ITだったり、仕事で覚えた専門知識の備忘録