FM ステレオ復調(Julia)

## 概要
ソフトウェア無線によるFM受信の技術メモです。今回は懸案であったコンポジットステレオ信号のソフトウェア復調をおこなってみました。
関連記事）
[ADALM PLUTOでFM受信してみる](https://qiita.com/0x20FE/items/0e229a4c31889903d019)(python版)
[ADALM PLUTOでFM受信(Julia版)](https://qiita.com/0x20FE/items/9b6498083ed52954aa42)

## 解説
コンポジット信号はL＋R信号、１９ｋHzのパイロット信号、パイロット信号に位相同期した38kHzでDBS変調されたL-R信号、さらにDARCの信号が多重されています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/115334/341ae653-8643-120c-c376-56f8ed5e6c87.png)
図　コンポジット信号のスペクトル例

復号の手順は、受信サービスさんに連載されている、「shu-chanの放送ネットワーク道しるべ」の手順に従って実装しました。

（参考）Shu-chan,[＜ＦＭ放送用受信機（その３）＞](https://www.jushin-s.co.jp/michi/download/88_n21.pdf),放送ネットワーク道しるべ中山道（小田井宿）,受信サービス株式会社

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/115334/6e5d9fe4-06ff-0784-c1ed-4e2b236a770c.png)

具体的手順は次の通り。
1. まずパイロット信号を再生
2. パイロット信号を低減
2. Fc=15kHz程度のLPFでL＋R信号をとりだす
2. 3KHzの信号をＭＩＸしたのちLPFをかけてL-R信号をとりだす。
2. L+R,L-RからL,Rを取り出し、De-empasisフィルタをかける


### パラメータとレジスタ格納領域を作成
サンプリングパラメータ、NCO設定データ、PLLパラメータとレジスタ、低域通過フィルタのレジスタ、ディエンファシスパラメータの格納領域を定義する。

```julia:class定義
mutable struct StereoDecoder
    FS_AUD # sampling rate of audio signal
    FS_DEM # sampling rate of composit signal(FS_AUD x 4)
    rNco::Nco # NCOクラス

    k1_pll # PLL parameter
    k2_pll
    z_pll  # PLL register
    
    aStatCicSin  # registers for CIC LPF
    aStatCicCos  
    aStatCicP
    aStatCicM
    
    a_deemp # ディエンファシス　パラメータ
    b_deemp
end
```
### パラメータ設定
NCOをパイロット信号の周波数19kHzで初期化し、ディエンファシスを50μ秒とする。
（参考）ディエンファシスフィルタの設計は、[ADALM PLUTOでFM受信してみる](https://qiita.com/0x20FE/items/0e229a4c31889903d019)(python版)

```julia:パラメータ設定
function cStereoDecoder_create(FsAud,FsDem)::StereoDecoder
    StereoDecoder(FsAud, FsDem, cNco_create(FsDem,19000.0),
                    5,2,0.0,
                    [],[],[],[],
                    [FsAud/10000+1,-(FsAud/10000-1)],[1,1])
end
```
### デコード処理
処理手順は
1. NCOを19kHzのパイロット信号に同期させる
このため、NCOの正弦波出力とコンポジット信号をMIXしDC成分をとりだし、これが０になるようにNCOを制御する。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/115334/08c71973-40d0-bad8-fe3a-9aefeb76f175.png)

1. パイロット信号の振幅を計算する
NCOから余弦波を作成し、コンポジット信号にMIXしDC成分を取り出すと、パイロット信号の振幅の1/2が得られる
1. コンポジット信号のパイロット信号を打ち消す
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/115334/118a6d4a-18ae-1bca-17b9-ecc8026c88da.png)
※あまりきえていませんね。

1. 前記信号にLPFをかけて15kHzまでのL+R信号をとりだす
1. 前記信号にパイロット信号に同期した38kHzの信号をMIXし、LPFをかけてL-R信号をとりだす。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/115334/115a79d6-d298-0990-82bb-c8449bb1321a.png)

1. L＋R,L-RからL,R信号をとりだす
1. L,Rそれぞれにディエンファシスをかける
 ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/115334/7647a97d-6440-26ab-bdd8-9477c76983f3.png)


```julia:
function mStereoDecoder_proc!(rSD::StereoDecoder,x,nDiv)
    nSmp=size(dem)[1]
    nDur=nSmp ÷ nDiv
    nR=128;nCasc=4
    adm2=[]; amx2=[]    # storages for output
    for iPos=1:nDur:nSmp
        #NCOで位相を計算してから、19kHzのSin信号を作成する
        aPh = mNco_generatePhase!(rSD.rNco,nDur)
        aSin= cNco_sinWave(aPh)

        ### LOCK PILOT ###
        #--- mix ---
        mx= x[iPos:iPos+nDur-1] .* aSin
        #--- LPF ---
        (outSin,rSD.aStatCicSin)=cic_ma(mx,nR,nCasc,rSD.aStatCicSin)
        s= sum(outSin)/nDur *10 # 1/10%  *rSD.FS_DEM /2/75000 *10
        #--- PLL ---
        fctr = rSD.k1_pll * s + rSD.z_pll
        rSD.z_pll += rSD.k2_pll*s
        rSD.rNco.dFreqDelta =  -fctr

        ### CANCEL PILOT ###
        aCos= cNco_cosWave(aPh)
        mx= x[iPos:iPos+nDur-1] .* aCos
        #--- LPF ---
        (outCos,rSD.aStatCicCos)=cic_ma(mx,nR,nCasc,rSD.aStatCicCos)
        #--- amplitude of Pilot
        p= sum(outCos)/nDur  *2 #   *rSD.FS_DEM /2/75000 *10

        #--- cancel pilot form dem
        dm2= x[iPos:iPos+nDur-1] - aCos * p

        ### demod L-R ###
        aSin2=cNco_sin2Wave(aPh)

        mx2 = dm2 .* aSin2
        if size(amx2)[1]==0
            adm2=copy(dm2)
            amx2=copy(mx2)
        else
            append!(adm2,dm2)
            append!(amx2,mx2)
        end
    end
    
    (outP,rSD.aStatCicP)=cic_ma(adm2,4,4,rSD.aStatCicP)
    (outM,rSD.aStatCicM)=cic_ma(amx2,4,4,rSD.aStatCicM)
    outPd=outP[1:4:end]
    outMd=outM[1:4:end]
    Lr=outPd+outMd
    Rr=outPd-outMd
    L=filt(rSD.b_deemp,rSD.a_deemp,Lr)
    R=filt(rSD.b_deemp,rSD.a_deemp,Rr)
    return (L,R)
end
```

## 性能、残課題など
SNがよいFM局であれば、ほぼ問題なくステレオ分離できているが、信号が弱いとノイジーになります。このあたりは、PLLのパラメータ、分割区間の調整で改善するとよいのですが。

L＋R,L-R信号を取り出すためのダウンサンプル用LPFをCIC_Filterで代用しています。このため通過帯域がだれています。本線なので、FIR構成のフィルタなどできちんと処理した方がよいかと。

NCO,CICフィルタの実装解説は、またの機会に。