【python】各種変調方式(QPSK, 16QAM, 64QAM)でAWGN通信を行ったときのBERの比較

過去の自分の記事をみると、MATLABで全く同じことをしていましたが、もしかしたらpythonのほうが需要があったりするのではないかと思い作った次第です。

本記事の内容としては。

https://qiita.com/Seiji_Tanaka/items/a059090ee020a5d5c648
に書いたビット→シンボルの変調用関数が入ったpythonファイル「QAM_modulation」と、

https://qiita.com/Seiji_Tanaka/items/b29d8a5a3447680370ce
に書いたシンボル→ビットの復調用関数が入ったpythonファイル「QAM_demodulation」と

import numpy as np

def AWGN_noise(EbNo, N, sym_num, bit_num):
    """
    EbNo : 発生させたい雑音の強度を決める
    N : 受信機側アンテナの数
    sym_num : 送信シンボル系列長
    bit_num : 送信ビット系列長
    """

    #変調多値数を求める　QSPK : 2, 16QAM : 4, 64QAM : 6となる値
    QAM_order = int(bit_num / sym_num)
    SNR = QAM_order * 10 ** (EbNo / 10)
    No = 1 / SNR

    noise = np.random.normal(0, np.sqrt(No / 2), (N, sym_num)) \
            + 1j * np.random.normal(0, np.sqrt(No / 2), (N, sym_num))

    return noise

という3つのpythonファイルをインポートして、通信シミュレーションを行います。

import numpy as np
import matplotlib.pyplot as plt

# https://qiita.com/Seiji_Tanaka/items/a059090ee020a5d5c648
import QAM_modulation

# https://qiita.com/Seiji_Tanaka/items/b29d8a5a3447680370ce
import QAM_demodulation

# Qiitaのこのページの上部に記載
import channel_generation

# 送信側アンテナ数
M = 1
# 受信側アンテナ数
N = 1
# 送信ビット数
bit_num = 12 * 10 ** 5

# 計算するEbNoの範囲　freqは計測点の間隔の狭さに対応する
EbNo_min = 0
EbNo_max = 20
freq = 21
power_range = np.linspace(EbNo_min, EbNo_max, freq)

BER_QPSK = np.zeros((1, freq))
BER_16QAM = np.zeros((1, freq))
BER_64QAM = np.zeros((1, freq))
count = 0

for EbNo in power_range:
    TX_bit = np.random.randint(0, 2, (M, bit_num))

    # 変調
    TX_QPSK = QAM_modulation.bit_to_QPSK(TX_bit)
    TX_16QAM = QAM_modulation.bit_to_16QAM(TX_bit)
    TX_64QAM = QAM_modulation.bit_to_64QAM(TX_bit)

    # AWGN雑音を付加
    noise = channel_generation.AWGN_noise(EbNo, N, len(TX_QPSK[0, :]), bit_num)
    RX_QPSK = TX_QPSK + noise

    noise = channel_generation.AWGN_noise(EbNo, N, len(TX_16QAM[0, :]), bit_num)
    RX_16QAM = TX_16QAM + noise

    noise = channel_generation.AWGN_noise(EbNo, N, len(TX_64QAM[0, :]), bit_num)
    RX_64QAM = TX_64QAM + noise

    # 判定・復調
    RX_bit_QPSK = QAM_demodulation.QPSK_to_bit(RX_QPSK)
    RX_bit_16QAM = QAM_demodulation.QAM16_to_bit(RX_16QAM)
    RX_bit_64QAM = QAM_demodulation.QAM64_to_bit(RX_64QAM)

    # BERの計算
    error_sum_QPSK = np.sum(np.abs(TX_bit - RX_bit_QPSK))
    error_sum_16QAM = np.sum(np.abs(TX_bit - RX_bit_16QAM))
    error_sum_64QAM = np.sum(np.abs(TX_bit - RX_bit_64QAM))

    BER_QPSK[0, count] = error_sum_QPSK / (M * bit_num)
    BER_16QAM[0, count] = error_sum_16QAM / (M * bit_num)
    BER_64QAM[0, count] = error_sum_64QAM / (M * bit_num)
    count += 1

plt.semilogy(power_range, BER_QPSK[0,:], label="QPSK")
plt.semilogy(power_range, BER_16QAM[0,:], label='16QAM')
plt.semilogy(power_range, BER_64QAM[0,:], label="64QAM")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=0, fontsize=18)
plt.ylim(10**(-5), 10**(-0))
plt.xlabel('Eb/N0　[dB]')
plt.ylabel('BER')
ax = plt.gca()
ax.set_yscale('log')
plt.grid(which="both")
plt.show()

出力結果は、

という感じで得られると思います。

もしpythonで通信系の理論シミュレーションを行おうとする方で変復調用の関数を作るのに手間取っているひとがいたらぜひ使ってください。
こういうことで悩む時間というのは本当にもったいないと思いますので。