ADALM PLUTOでFM受信してみる

## 概要
ソフトウエア無線によるFM受信の技術メモです。RFアンプとAD変換部分は、Analog DevicesのPLUTOを使い、Windows環境のpythonにてFM受信を行ってみました。

※更新
 - 遅延検波の解説を更新(2021/9/23)
 - ディエンファシスの追加(2022/6/19)　＜-忘れてた。
 - 読みにくかったのでリライトしました(2022/8/12,2023/9/14)
- [NEW] 実行環境の準備に、PLUTOの環境変数の設定を追加(2024/6/28) 
* 関連記事
[ADALM-PLUTOを有線LANに接続する](https://qiita.com/0x20FE/items/8e5daeadeed327f9737e) (2022/8/12)

#### 想定するアプリケーションと実行環境
想定するアプリケーション：SDR、FM受信、音声処理
実行環境：python, windows10, ADALM PLUTO
他のデバイスでも受信処理とオーディオ出力の部分はそのまま使えるかと。

#### 実行環境の準備
PLUTOを使うために次のリンクに従ってWindows版のLibiio,python,pyadi-iioをインストールする。

https://wiki.analog.com/resources/tools-software/linux-software/pyadi-iio

初期設定では、ローカル周波数の設定範囲が(325 - 3800MHz）であるため、（teratermなど使って）シリアル通信でログインし、環境変数を設定しrebootすることで、周波数設定範囲を(70 - 6000 MHz)に拡大する。
（参考）[Customizing the Pluto configuration](https://wiki.analog.com/university/tools/pluto/users/customizing) Updating to the AD9364
```:シリアルコンソールから
# fw_setenv attr_name compatible
# fw_setenv attr_val ad9364
# reboot
```

## 解説
### FM検波理論（遅延検波）
瞬時周波数成分を取り出す（FM検波、FM復調）方法として遅延検波を利用します。

s(t)を伝送信号としたとき、ベースバンドFM信号$x(t)$の瞬時周波数f(t)は
$f(t) = K s(t)$ Kは定数
ベースバンドFM信号$x(t)$をサンプリング間隔Tでサンプルすると
$x(t) = exp(j2πft)$
$x_n=exp(j2πnfT)$（fは瞬時周波数）

この信号を１サンプル遅延させ複素共役した値と積算する。
$ g_n = x_n  x_{n-1}^* $
$       = exp(j2πnfT)exp(-j2πnfT+j2πfT)$
$       = exp(j2πfT)$

したがって、np.angle()で複素数の偏角を求め$2πT$で割ると瞬時周波数fが得られる。
$f = \frac{1}{2πT} np.angle(x_n x_{n-1}^*)$
       

```python:遅延検波
# (i) bbIQ
# (o) dem
dem = np.angle(bbIQ[1:] * np.conj(bbIQ[:-1])) / math.pi 
```

### サンプル間隔の選定
サンプル間隔$T$が大きくなると前記$g_n$の偏角が$±π$を超え、計算結果が不連続となり異常値が発生する。逆にサンプル間隔$T$が小さいすぎると偏角の値が小さくなり、量子化雑音の影響が大きくなる。このためサンプル間隔Tは注意深く選定する必要がある。
$ 2πTf_{max} < 2π$ となる$T$は、 $ T < \frac{1}{ f_{max}} $

### ディエンファシスフィルタの導出
日本のFM放送はτ=50usのエンファシスがかかっているため、ディエンファシスフィルタをかける。
$He(s) = 1 + τs $ エンファシス 
$Hd(s) = \frac{1}{1 + τs}$ディエンファシス


デジタルフィルタ構成とするためHd(s)に双一次変換を施す。
$s = \frac{2}{T_s } \frac{1 - z^{-1}}{ 1 + z^{-1}}$ 

サンプリング 周波数を44100Hzとすると
$Hd(z) = \frac{ 1 - z^{-1} }{5.41 - 3.4z^{-1} }$ 

```python:ディエンファシスフィルタ
# (i) aud
# (o) demp
deemp_n = np.array([1,1])*2
deemp_d = np.array([5.41,3.41])
demp = scipy.signal.lfilter(deemp_n,deemp_d,aud)
```

### Pythonでのサウンド出力（PyAudio）
別記事を参照してください。

https://qiita.com/0x20FE/items/d8a4a477855b6e3aa36f

```python:pyaudio
#初期化
pyAud = pyaudio.PyAudio()
#サンプリング周波数、チャンネル、フォーマット
sndStrm = pyAud.open(rate=SND_FS,channels=SND_CH_OUT,format=pyaudio.paInt16,input=False,output=True)

# 途中の繰り返し処理
# (i) pack
# (o) sndStrm
# (g) SND_CHUNK バッファ長
while pack.size > SND_CHUNK :
	# ±1に正規化された信号をSND_CHUNK毎に16bit固定少数変換して出力
    aBuf = pack[0:SND_CHUNK]*8000
    sndBuf = aBuf.astype(np.int16).tobytes()
    sndStrm.write(sndBuf)
    pack = pack[SND_CHUNK:]
    
    
#後始末    
sndStrm.stop_stream()
sndStrm.close()
pyAud.terminate()
```

### サンプルプログラム

1. sdr=SdrInit()にてPlutoの設定を行い、オブジェクトを得る。
2. sdr.rx()でサンプルを取り出し、demLoop()にて遅延検波処理を行い音声データを得る。
3. 遅延検波の後で、音声デバイスのサンプルレートにダウンサンプルし、デエンファシスをかける。
4. メインループでは、demLoopで取り出した音声データを、一旦バッファ(pack)に格納し、pyaudioのチャンクバッファサイズよりデータが溜まったところで、チャンクサイズごとに取り出し、pyaudio用のバイナリデータに変換してpyaudioのデバイスに（ブロックモードで）書き込む。

```python:PlutoFM.py
import adi                          #pluto 
import matplotlib.pyplot as plt     #graph
import numpy as np
import scipy.signal                 # US/DS/Filter
import math                         # pi.sic.cos
import pyaudio

# System parameters
################################################
#-- application
################################################
SND_FS = 44100 # 48000
SND_CHUNK = 1024*2
SND_CH_OUT = 1

sample_rate_audio = SND_FS
################################################
#-- baseband IQ parameter
################################################
#(i) sample_rate_audio
decimate_demod = 5 
sample_rate_demod = sample_rate_audio * decimate_demod
################################################
#-- radio interface param
################################################
# PLUTO
# LO 325MHz to  3800MHz
# BW < 20MHz
# FS 521ksps to 61.4Msps # cat spec is 65.2Ksps to  61.4Msps
################################################
#(i) sample_rate_demod
decimate_capture = 4 
sdrSample_rate = sample_rate_demod * decimate_capture

sdrBufSize = 1024*20
sdrFcenter = 82500000 # 中心周波数82.5MHzの場合
#sdrFcenter = 80000000
sdrBandWidth = 400000
################################################
# Create radio object
def SdrInit():
    # (i) 
    # (o) radio Object
    sdr = adi.Pluto()
    sdr.rx_rf_bandwidth = sdrBandWidth
    sdr.rx_lo = sdrFcenter
    sdr.sample_rate = sdrSample_rate
    sdr.rx_buffer_size=sdrBufSize
    #sdr.tx_lo = 2000000000
    #sdr.tx_cyclic_buffer = True
    #sdr.tx_hardwaregain = -30
    sdr.gain_control_mode = "slow_attack"
    # Read back properties from hardware
    #print(sdr.rx_hardwaregain)
    # Get complex data back
    sdr.rx_enabled_channels = [0]
    return sdr

#=========================================================
#　FM　復調
#=========================================================
deemp_n = np.array([1,1])*2
deemp_d = np.array([5.41,3.41])
demLastIQ=0
def demLoop(demLastIQ):
    smp = sdr.rx()
    decim = scipy.signal.decimate(smp,decimate_capture)
    bbIQ = np.insert(decim,0,demLastIQ)
    #bbIQ = decim
    demLastIQ = decim[-1]
    #------------------------
    # 遅延検波
    # -+-[Delay]-- [x]---[ angle ]-->
    #  |            ^
    #  +-[Conj ]----+
    #------------------------
    dem = np.angle(bbIQ[1:] * np.conj(bbIQ[:-1])) / math.pi 
    #  decimate to audio sample rate
    aud = scipy.signal.decimate(dem,decimate_demod)
    #deEmphasis
    demp = scipy.signal.lfilter(deemp_n,deemp_d,aud)
    return demp,demLastIQ

if __name__ == "__main__":
    sdr = SdrInit()
    #====================================================
    #--- setup Audio
    #====================================================
    pyAud = pyaudio.PyAudio()
    sndStrm = pyAud.open(rate=SND_FS,channels=SND_CH_OUT,format=pyaudio.paInt16,input=False,output=True)
    #====================================================
    #--- setup FIFO
    #====================================================
    pack = np.array([])
    #====================================================
    #--- FM受信ループ
    #====================================================
    for cnt in range(2000): # while: でもよい
        #====================================================
        #--- チャンクサイズ以上たまったらサウンドデバイスに出力
        #====================================================
        while pack.size > SND_CHUNK :
            aBuf = pack[0:SND_CHUNK]*8000
            sndBuf = aBuf.astype(np.int16).tobytes()
            sndStrm.write(sndBuf)
            pack = pack[SND_CHUNK:]
            #plt.plot(aBuf,linestyle='None',marker="x")
            #plt.psd(aBuf,NFFT=1024,Fs=sample_rate_audio , Fc=0)
        #====================================================
        #--- 新たなIF信号を受信し音声データに復調しFIFOにつめる
        #====================================================
        aud,demLastIQ = demLoop(demLastIQ)
        pack = np.append(pack, aud)
    #====================================================
    # teardown 処理
    #====================================================
    sndStrm.stop_stream()
    sndStrm.close()
    pyAud.terminate()

```
## 結言
パラメータ調整は必要ですが、遅延検波で簡単にFMを復調できます。