PythonでFM信号の直交復調を試す

pythonに慣れるのと、SDRにおいてIQ復調(Quadrature demodulation)の理解をする。

参考
直交復調の解説
http://einstlab.web.fc2.com/IQ/IQ.pdf

FM変調波を用意する

サンプリング周波数は1MHz
キャリア周波数は20kHz
信号周波数は1kHz
ナイキスト周波数はサンプリング周波数の半分

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

sample_rate = 1.0e6
carrier_freq = 20e3
signal_freq = 1e3
fn = sample_rate/2.0
duration = 1.0/signal_freq
xs = np.arange(0, duration, 1/sample_rate)
ts = xs * 1000

変調するまで

# signal
ys1 = np.sin(2.0*np.pi*signal_freq*xs)
plt_graph(plt, "time[ms]", "Signal", ts, ys1)

# carrier
carrier = np.sin(2.0*np.pi*carrier_freq*xs)
plt_graph(plt, "time[ms]", "Carrier Singal", ts, carrier)

# modulation
# m = 3.5
m = 1
ys2 = np.sin(2.0*np.pi*carrier_freq*xs + m*ys1)
plt_graph(plt, "time[ms]", "FM Singal", ts, ys2)

変調指数mはとりあえず1にしておく。

変調指数が小さいと変調信号はキャリア信号と違いはそこまでわからない。

直交復調をする

IQ信号を生成するためのLOを計算する。(いわゆるNumerical Controlled Oscillatorのようなもの)

# LO generation
lo_i = np.cos(2.0*np.pi*carrier_freq*xs)
lo_q = np.sin(2.0*np.pi*carrier_freq*xs)
plt.figure()
plt.xlabel("time[ms]")
plt.title("LO Signal")
plt.plot(ts, lo_i, 'r')
plt.plot(ts, lo_q, 'g')

ミキシング

# mixing
bbi = ys2*lo_i
bbq = ys2*lo_q
plt.figure()
plt.xlabel("time[ms]")
plt.title("Mixed Signal")
plt.plot(ts, bbi, 'r')
plt.plot(ts, bbq, 'g')

LPFを通す

実際のダイレクトコンバージョンではデシメーション後にLPFを通す。

# LPF
fp = 1.1e3
# fs = 2e6
# gpass = 1
# gstop = 40

Wp = fp/fn
# Ws = fs/fn
fir = signal.firwin(127, Wp)
bbi2 = signal.lfilter(fir, 1, bbi)
bbq2 = signal.lfilter(fir, 1, bbq)

plt.figure()
plt.xlabel("time[ms]")
plt.title("LPF Signal")
plt.plot(ts, bbi2, 'r')
plt.plot(ts, bbq2, 'g')

LPFのタップ数は適当

参考にしているPDF p13に近い形がでてくる。
(p12のベースバンド波形からはLPFを通しただけではp13の波形は出てこないが...)

位相を計算する

# arctan
phi_sig = np.ndarray(0)

for index, (i, q) in enumerate(zip(bbi2, bbq2)):
    phi_sig = np.append(phi_sig, np.arctan2(q, i))

plt.figure()
plt.xlabel("time[ms]")
plt.title("phi")
plt.plot(ts, phi_sig, 'g')

諸事情あってforで計算している

差を取る

# diff
demod = np.diff(phi_sig)
plt_graph(plt, "time[ms]", "Demodulated", ts[:len(demod)], demod)

LPFの計算で最初のところが連続ではないため、位相がおかしいので長く計算して真ん中を見る。
durationを3倍にする。

start = int(len(xs)/3)
print(start)
end = start + int(len(xs)/3)
print(end)
# plt_graph(plt, "time[ms]", "Demodulated", ts[:len(demod)], demod)
plt_graph(plt, "time[ms]", "Demodulated", ts[start:end], demod[start:end])

cosと言えなくもない波が出てきたが、波形が小さいし、どうにも高周波がのってしまっている。
微分の仕方がよくないのかもしれない。

変調指数を変えてみる

mを5ぐらいにすると、PDFのp12のような波がでてくる。

しかし、これをLPFに通したものは、p13とだいぶ違う上に、arctan2で折り返してしまう。

この折り返しをIとQの符号をみて何となく修正した。

th_q = list(map(lambda x: 1 if x > 0 else 0, bbq2))
th_i = list(map(lambda x: 1 if x > 0 else 0, bbi2))

# arctan
phi_sig = np.ndarray(0)
hoge = np.diff(th_q)
hoge = np.append(hoge, 0)
offset = 0
for index, (i, q) in enumerate(zip(bbi2, bbq2)):
    phi_sig = np.append(phi_sig, np.arctan2(q, i) + offset)

    if th_i[index] == 0 and index > 0:
        if hoge[index] > 0:
            offset = offset - np.pi * 2.0
        elif hoge[index] < 0:
            offset = offset + np.pi * 2.0