【調變解調】FM 調頻

說明

學習數位訊號處理演演算法時整理的學習筆記。同系列文章目錄可見 《DSP 學習之路》目錄，程式碼已上傳到 Github - ModulationAndDemodulation。本篇介紹 FM 調頻訊號的調變與解調，內附全套 MATLAB 程式碼。

1. FM 調變演演算法

1.1 FM 訊號描述

用調變訊號去控制載波的瞬時頻率，使其按照調變訊號的規律變化，當調變訊號是模擬訊號時，這個過程就被稱為調頻（FM）。FM 訊號的時域表示式為：

式中：\(A\) 為載波恆定振幅，\(K_f\) 為調頻靈敏度（單位 \(rad/(s{\cdot}V)\)），\(m(t)\) 是調變訊號（攜帶要發出去的資訊），\(cos{\omega_ct}\) 是載波，\(\omega_c\) 是載波角頻率，與載波頻率 \(f_c\) 之間的關係為 \(\omega_c=2{\pi}f_c\)。由式 \((1)\) 可得 FM 訊號相對於載頻 \({\omega_c}\) 的瞬時頻偏為：

由 \((1)\) 式可知 FM 訊號相對於載波相位 \({\omega}_ct\) 的瞬時相位偏移隨 \(m(t)\) 的積分呈線性變化，由 \((2)\) 式可知 FM 訊號相對於載頻 \({\omega_c}\) 的瞬時頻率偏移隨 \(m(t)\) 成線性變化，比例係數都為 \(K_f\)。有時候會用 \(k_f\) 表示調頻靈敏度（單位 \(Hz/V\)），它與 \(K_f\) 之間的關係為 \({K_f}=2{\pi}{k_f}\)，需注意這個轉換關係。FM 的調頻指數（調變指數） \({\beta}\) 定義如下，其中 \(W\) 是基頻訊號（調變訊號）\(m(t)\) 的頻寬（或最高頻率）：

若 \(m(t)\) 為單一頻率的正弦波（即 \(m(t)={A_m}cos({2{\pi}{f_m}t})\) 時，此時 \(W={f_m}\)），則調變指數的表示式如下，調變指數的值正好與最大相位偏移相同，其中 \({\Delta}f={\beta}{f_m}\) 表示最大頻偏。

FM 訊號的頻域表示式比較複雜，下面分成窄頻調頻與寬頻調頻兩種情況來簡化討論。

（1）窄頻調頻

一般將由 \(m(t)\) 引起的最大瞬時相位偏移遠小於 \(30^{\circ}\) 的情況稱為窄頻 FM，即滿足以下條件時，FM 訊號的頻譜寬度比較窄，稱為窄頻調頻（NBFM）。窄頻調頻佔據的頻寬較窄，傳輸資料量有限，主要應用於無線語音的傳輸（比如無線對講機）。

此時，式 \((1)\) 可以近似為：

對其做傅立葉變換，得到窄頻調頻（NBFM）訊號的頻譜（幅度譜）表示式：

式中，\(M(\omega)\) 是調變訊號 \(m(t)\) 的頻譜。與 AM 訊號不同的是，NBFM 訊號的兩個邊頻分別乘了因式 \(1/(\omega-\omega_c)\) 和 \(1/(\omega+\omega_c)\)，由於因式是頻率的函數，所以這種加權是頻率加權，加權的結果引起調變訊號頻譜的失真，此外， NBFM 的一個邊帶和 AM 反相。

（2）寬頻調頻

當不滿足 \((5)\) 式所表示的條件時，FM 訊號被稱為寬頻調頻（WBFM）。寬頻調頻所佔用的頻頻寬度比較寬，傳輸資料量大，主要應用於調頻立體聲廣播。WBFM 的時域表示式無法進一步簡化，首先考慮單音調變的情況，然後把分析的結論推廣到多音調變，當 \(m(t)={A_m}cos({{\omega_m}t})\) 時，帶入式 \((1)\) 展開可得：

式中 \(J_n(\beta)\) 為第一類 \(n\) 階貝塞爾函數，它是調頻指數 \(\beta\) 的函數。對其做傅立葉變換，得到單音調變時寬頻調頻（WBFM）訊號的頻譜（幅度譜）表示式：

由式 \((8)\) 和式 \((9)\) 可見，WBFM 調頻訊號的頻譜由載波分量 \(\omega_c\) 和無數邊頻 \({\omega_c}\pm{n\omega_m}\) 組成。當 \(n=0\) 時是載波分量 \(\omega_c\)，其幅度為 \(AJ_0(\beta)\)，當 \(n\not=0\) 時就是對稱分佈在載頻兩側的邊頻分量 \({\omega_c}\pm{n\omega_m}\)，其幅度為 \(AJ_n(\beta)\)，相鄰邊頻之間的間隔為 \(\omega_m\)。且當 \(n\) 為奇數時，上下邊頻極性相反，當 \(n\) 為偶數時極性相同。由此可見，FM 訊號的頻譜不再是調變訊號頻譜的線性搬移，而是一種非線性過程。

1.2 FM 訊號的頻寬與功率分配

寬頻調頻訊號的頻譜包含無窮多個頻率分量，因此理論上調頻訊號的頻頻寬度為無限寬。但是，實際上邊頻幅度 \(J_n(\beta)\) 隨著 \(n\) 的增大而逐漸減小，因此只要取適當的 \(n\) 值使邊頻分量小到可以忽略的程度，調頻訊號可以近似認為具有有限頻譜。通常採用的原則是：訊號的頻頻寬度應包括幅度大於未調載波的 \(10\%\) 以上的邊頻分量，即 \(\lvert{J_n(\beta)}\rvert\geq0.1\)。根據經驗，當 \(\beta\geq1\) 時，取邊頻數 \(n=\beta+1\) 即可，因為 \(n>\beta+1\) 以上的邊頻幅度 \(J_n(\beta)\) 均小於 \(0.1\)，相應產生的功率均在總功率的 \(2\%\) 以下，可以忽略不計，根據這條經驗規則，調頻波的有效頻寬為：

這就是廣泛用於計算調頻訊號頻寬的卡森（Carson）公式。式中 \(\beta\) 為調頻指數，\(W\) 為基頻訊號（調變訊號）\(m(t)\) 的頻寬（或最高頻率）。對於單音調變訊號（即 \(m(t)={A_m}cos({2{\pi}{f_m}t})\) 時，此時 \(W={f_m}\)），將 \((4)\) 式帶入 \((10)\) 式，可得：

• 當 \(\beta \ll 1\) 時（對應窄頻調頻），頻寬可近似估計為 \(B_{NBFM}\approx2f_m\)，這時，頻寬由第一對邊頻分量決定，頻寬只隨調變頻率 \(f_m\) 變化，而與最大頻偏 \({\Delta}f\) 無關。
• 當 \(\beta \gg 1\) 時（對應寬頻調頻），頻寬可近似估計為 \(B_{WBFM}\approx2{\Delta}f\)，這時，頻寬由最大頻偏 \({\Delta}f\) 決定，而與調變頻率 \(f_m\) 無關。

利用帕塞瓦爾定理以及貝塞爾函數性質，不難求得 FM 訊號功率 \(P_{FM}\) 與未調載波功率 \(P_c\) 的關係如下：

上式說明，調頻訊號的平均功率 \(P_{FM}\) 等於未調載波的平均功率 \(P_c\)，即調變後總的功率不變，只是將原來載波功率中的一部分分配給每個邊頻分量，所以，調變過程只是進行功率的重新分配，而分配的原則與調頻指數 \(\beta\) 有關。

1.3 FM 訊號的調變方法

FM 訊號的調變方法有 3 種，一種是直接調頻法，一種是間接調頻法，第三種是正交調變法。

（1）直接調頻法

直接調頻法是用調變訊號 \(m(t)\) 直接控制高頻振盪器，讓迴路元件的引數發生改變，使其輸出頻率按調變訊號的規律線性地變化，常用的元件是變容二極體。直接調頻法的主要優點是在實現線性調頻的要求下，可以獲得較大的頻偏，且實現電路簡單；主要缺點是頻率穩定度不高，往往需要採用自動頻率控制系統來穩定中心頻率。

（2）間接調頻法

間接調頻法也被稱為倍頻法、阿姆斯特朗法。該方法先將調變訊號積分，然後對載波進行調相，即可產生一個 NBFM 訊號，再經 \(n\) 次倍頻器得到 WBFM 訊號，間接調頻法的優點是頻率穩定性好，缺點是需要多次倍頻和混頻，電路較複雜。如下圖所示，根據式 \((6)\) 窄頻調頻的時域表示式可知，虛線框內所示部分可以用來獲得 NBFM 訊號。

上圖中的倍頻器是一個非線性器件，以理想平方律器件為例，其輸出 \(x_o(t)\) 與輸入 \(x_i(t)\) 的關係為 \({x_o}(t)=a{x_i^2}(t)\)，其中 \(a\) 為常數，將式 \((6)\) 帶入，得到：

接著將 \(x_o(t)\) 經過一個帶通濾波器，就可以將其中的直流分量濾除，從而得到一個新的 FM 訊號 \(x'_o(t)\)，如下：

這個訊號的載頻和相位偏移均增為原來的 \(2\) 倍，由於相位偏移增為 \(2\) 倍，因而調頻指數也必然增為 \(2\) 倍，同理，經 \(n\) 次倍頻後可以使調頻訊號的載頻和調頻指數增為 \(n\) 倍。增大調變指數時，一般需要保持載頻不變或者不增大太多，這個時候可以接著使用一個混頻器配一個帶通濾波器來進行下變頻，混頻器只改變載頻而不影響頻偏，具體實現方法可參考《通訊原理》相關資料。

（3）正交調變法

將 \((1)\) 式進行三角展開，可以得到：

正交調變流程如下：

1. 對調變訊號 \(m(t)\) 進行積分，得到 \(\Phi=K_f\int_0^{t}{m(\tau)}d{\tau}\)。
2. 對積分後的訊號分別取餘弦和正弦，得到 \(I\) 路資料 \(I(t)=cos(\Phi)\) 與 \(Q\) 路資料 \(Q(t)=sin(\Phi)\)。
3. 分別乘以載波 \(Acos(\omega_ct)\) 與 \(-Asin(\omega_ct)\) 後相加，得到 FM 訊號 \(s_{FM}(t)=I(t)Acos(\omega_ct)-Q(t)Asin(\omega_ct)\)。

其中第三步也可以先將 \(I(t)\) 和 \(Q(t)\) 組成一個覆信號 \(Z(t)=I(t)+iQ(t)\)，然後乘以覆載波 \(exp(i\omega_ct)\)，最後取實部，即 \(s_{FM}(t)=Real\left[Z(t)exp(i\omega_ct)\right]\)，兩種方法是等價的。

1.4 窄頻 FM 訊號範例

上一小節中介紹的調變方法都與硬體實現相關，接下來使用 MATLAB 軟體來直接生成 NBFM 訊號。調變訊號 \(m(t)\) 可以是確知訊號，也可以是隨機訊號。當 \(m(t)\) 是確知訊號時，不妨假設 \(m(t)\) 的時域表示式如下，對應的基頻頻寬 \(W=f_m\)：

各調變引數取值：\(A=1\)，\(f_m=2500Hz\)，\({\beta}=10^{-6}\)，\(f_c=20000Hz\)。訊號取樣率 \(f_s=8{f_c}\)，模擬總時長為 \(2s\)。FM 調變效果如下圖所示（為了美觀，時域只顯示前 500 個點），調變訊號 \(m(t)\) 雙邊幅度譜有四根離散譜線（\({\pm}2500Hz\)、\({\pm}1250Hz\)），調變訊號 \(m(t)\) 積分後的雙邊幅度譜有五根離散譜線（\(0Hz\)、\({\pm}2500Hz\)、\({\pm}1250Hz\)），NBFM 調頻訊號 \(s(t)\) 雙邊幅度譜有十根離散譜線（\({\pm}22500Hz\)、\({\pm}21250Hz\)、\({\pm}20000Hz\)、\({\pm}18750Hz\)、\({\pm}17500Hz\)）。NBFM 調頻訊號 \(s(t)\) 的頻寬滿足公式 \(B_{NBFM}\approx2f_m=5000Hz\)，程式碼詳見附錄 main_modFM_example1.m 與 mod_fm.m，設定 beta = 1e-6 即可。

當 \(m(t)\) 是隨機訊號時，不妨假設基頻訊號頻寬為 \(W={f_H}=3000Hz\)，各調變引數取值：\(A=1\)，\({\beta}=10^{-6}\)，\(f_c=20000Hz\)。訊號取樣率 \(f_s=8{f_c}\)，模擬總時長為 \(2s\)。FM 調變效果如下圖所示（為了美觀，時域只顯示前 500 個點），調變訊號 \(m(t)\) 雙邊幅度譜中間譜峰的範圍約為 \(-3000Hz{\sim}3000Hz\)，調變訊號 \(m(t)\) 積分後的雙邊幅度譜中間譜峰的範圍依然約為 \(-3000Hz{\sim}3000Hz\)，NBFM 調頻訊號 \(s(t)\) 雙邊幅度譜有兩根離散譜線（\({\pm}20000Hz\)）及兩個譜峰（範圍約為 \(-23000Hz{\sim}-17000Hz\)、\(17000Hz{\sim}23000Hz\)）。NBFM 調頻訊號 \(s(t)\) 的頻寬滿足公式 \(B_{NBFM}\approx2f_H=6000Hz\)，程式碼詳見附錄 main_modFM_example2.m 與 mod_fm.m，設定 beta = 1e-6 即可。

1.5 寬頻 FM 訊號範例

當 \(m(t)\) 是確知訊號時，設 \(f_m=2500Hz\)， \(\beta=4\)，其他引數與 1.4 節中的確知訊號相同，FM 調變效果如下圖所示，頻寬約為 \(B_{FM}=2(\beta+1)f_m=25000Hz\)，程式碼詳見附錄 main_modFM_example1.m 與 mod_fm.m，設定 fm = 2500, beta = 4 即可。

.當 \(m(t)\) 是確知訊號時，設 \(f_m=1000Hz\)， \(\beta=15\)，其他引數與 1.4 節中的確知訊號相同，FM 調變效果如下圖所示，頻寬約為 \(B_{WBFM}\approx2{\Delta}f=2{\beta}f_m=30000Hz\)，程式碼詳見附錄 main_modFM_example1.m 與 mod_fm.m，設定 fm = 1000, beta = 15 即可。

2. FM 解調演演算法

解調是調變的逆過程，其作用是從接收的已調訊號中恢復原基頻訊號（即調變訊號）。FM 解調的方法也分為相干解調和非相干解調，普通的相干解調僅適用於 NBFM 訊號，正交相干解調與非相干解調對 NBFM 訊號和 WBFM 訊號均適用。下面分別用幾種不同方法對 1.4 與 1.5 節中確知訊號的 FM 調變結果進行解調。

2.1 非相干解調（鑑頻器）

微分器和包絡檢波器是一種最常見的鑑頻器。其中微分器的作用是把幅度恆定的調頻波 \(s_{FM}(t)\) 變成幅度和頻率都隨調變訊號 \(m(t)\) 變化的調幅調頻波 \(s_d(t)\)，即：

包絡檢波器則將其幅度變化檢出並濾去直流，即得解調輸出 \(m_o(t)={K_d}{K_f}m(t)\)，式中 \(K_d\) 為鑑頻器靈敏度（單位 \(V/(rad/s)\)）。FM 非相干解調（鑑頻器）一般有以下四個步驟，其中第一步是微分器，第二至第四步是包絡檢波器，與 AM 包絡檢波的流程一樣。

1. 第一步：求微分，得到 \(s_d(t)\)。
2. 第二步：全波整流（對 \(s(t)\) 取絕對值）或半波整流（將 \(s(t)\) 小於 \(0\) 的地方置零）。
3. 第三步：低通濾波器濾除高頻載波，濾除 \(2{\omega}_c\) 或 \({\omega}_c\)。
4. 第四步：去除直流分量（減去自身均值）。

對 1.4 節中的 FM 訊號，設定 \(f_m=2500Hz\)，\({\beta}=10^{-6}\)，訊雜比 \(SNR=150dB\)，非相干解調效果如下，最大頻偏 \({\Delta}f={\beta}f_m=0.0025Hz\) 過小，解調效果很差。

對 1.5 節中的 FM 訊號，設定 \(f_m=2500Hz\)，\({\beta}=4\)，訊雜比 \(SNR=50dB\)，非相干解調效果如下，解調後幅度放大係數 \(k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx7.61\)，使用這個係數放大解調訊號幅值，然後計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0920\)。

對 1.5 節中的 FM 訊號，設定 \(f_m=1000Hz\)，\({\beta}=15\)，訊雜比 \(SNR=50dB\)，非相干解調效果如下，解調後幅度放大係數 \(k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx5.07\)，使用這個係數放大解調訊號幅值，然後計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0634\)。

非相干解調（鑑頻器）的程式碼詳見 lpf_filter.m、demod_fm_method1.m 和 main_demodFM_example1.m。

2.2 非相干解調（鑑頻器 - 希爾伯特變換法）

鑑頻器中的包絡檢波器可以通過希爾伯特變換法來實現，解調時無需任何載頻資訊，這樣，FM 非相干解調（鑑頻器）可以分為以下三個步驟。

1. 第一步：求微分，得到 \(s_d(t)\)。
2. 第二步：計算 \(s_d(t)\) 的希爾伯特變換，得到一個覆信號（實部為 \(s_d(t)\)，虛部為其希爾伯特變換結果），對所得覆信號取模，即為 \(s_d(t)\) 的包絡。
3. 第三步：去除直流分量（減去自身均值）。

對 1.4 節中的 FM 訊號，設定 \(f_m=2500Hz\)，\({\beta}=10^{-6}\)，訊雜比 \(SNR=150dB\)，非相干解調效果如下，最大頻偏 \({\Delta}f={\beta}f_m=0.0025Hz\) 過小，解調效果很差。

對 1.5 節中的 FM 訊號，設定 \(f_m=2500Hz\)，\({\beta}=4\)，訊雜比 \(SNR=50dB\)，非相干解調效果如下，解調後幅度放大係數 \(k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx4.87\)，使用這個係數放大解調訊號幅值，然後計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0492\)。

對 1.5 節中的 FM 訊號，設定 \(f_m=1000Hz\)，\({\beta}=15\)，訊雜比 \(SNR=50dB\)，非相干解調效果如下，解調後幅度放大係數 \(k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx3.26\)，使用這個係數放大解調訊號幅值，然後計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0433\)。

非相干解調（鑑頻器 - 希爾伯特變換法）的程式碼詳見 demod_fm_method2.m 和 main_demodFM_example2.m。

2.3 相干解調

相干解調時，為了無失真地恢復原基頻訊號，接收端必須提供一個與調變載波嚴格同步（同頻同相）的本地載波（稱為相干載波，可使用鎖相環技術得到）。普通的相干解調僅適用於 NBFM 訊號，其解調框圖如下：

其中：

相干載波 \(c(t)=-sin({\omega_c}t)\)，則相乘器的輸出為：

經低通濾波器取出其低頻分量：

再經微分器，即得解調輸出：

可見，相干解調可以恢復原調變訊號，這種解調方法與線性調變中的相干解調一樣，要求本地載波與調變載波同步，否則將使解調訊號失真。NBFM 相干解調可以分為以下幾步：

1. 第一步：乘以相干載波（即乘以 \(-sin({\omega_c}t+\phi_0)\)），注意載波初始相位。
2. 第二步：低通濾波濾除高頻載波，濾除 \(2\omega_c\)。
3. 第三步：求微分。

對 1.4 節中的 FM 訊號，設定 \(f_m=2500Hz\)，\({\beta}=10^{-6}\)，訊雜比 \(SNR=150dB\)，相干解調效果如下，解調後幅度放大係數 \(k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx222.30\)，使用這個係數放大解調訊號幅值，然後計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.1313\)。更改相干載波的初始相位為 \({\phi_0}=\pi/4,\pi/2\)，或者更改相干載波的中心頻率為 \(0.8f_c,1.2f_c\) 後，解調效果變差，說明這種方法對相干載波同頻同相的要求較高，魯棒性不夠強悍，可使用鎖相環技術來改善這一缺點。

對 1.5 節中的 FM 訊號，設定 \(f_m=2500Hz\)，\({\beta}=4\)，訊雜比 \(SNR=50dB\)，相干解調效果如下，可知這種方法不適用於非窄頻調頻訊號。

對 1.5 節中的 FM 訊號，設定 \(f_m=1000Hz\)，\({\beta}=15\)，訊雜比 \(SNR=50dB\)，相干解調效果如下，可知這種方法不適用於非窄頻調頻訊號。

相干解調的程式碼詳見 lpf_filter.m、demod_fm_method3.m 和 main_demodFM_example3.m。

2.4 數位正交解調

數位正交解調也屬於相干解調的一種，但這種方法具有較強的抗載頻失配能力，不要求相干載波嚴格的同頻同相。FM 數位正交解調一般有以下四個步驟：

1. 第一步：乘以正交相干載波得到 \({s_I}(t)\) 與 \({s_Q}(t)\)，即 \({s_I}(t)=s(t)cos({\omega_ct}+{\phi_0})\)，\({s_Q}(t)=-s(t)sin({\omega_ct}+{\phi_0})\)。
2. 第二步：低通濾波器濾除 \({s_I}(t)\) 與 \({s_Q}(t)\) 中的高頻分量。
3. 第三步：通過反正切計算相位 \(\Phi(t)=atan\left[\frac{s_Q(t)}{s_I(t)}\right]=k\int_0^{t}{m(\tau)}d{\tau}\)。
4. 第四步：對相位進行差分，得到解調結果 \(m_o(t)\)。

反正切運算在硬體中難以實現，通過一階近似，上面第三步與第四步可合併為以下式子：

對 1.4 節中的 FM 訊號，設定 \(f_m=2500Hz\)，\({\beta}=10^{-6}\)，訊雜比 \(SNR=150dB\)，數位正交解調效果如下，解調後幅度放大係數 \(k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx17782852.62\)，使用這個係數放大解調訊號幅值，然後計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.1318\)。

對 1.5 節中的 FM 訊號，設定 \(f_m=2500Hz\)，\({\beta}=4\)，訊雜比 \(SNR=50dB\)，數位正交解調效果如下，解調後幅度放大係數 \(k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx4.48\)，使用這個係數放大解調訊號幅值，然後計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0391\)。

對 1.5 節中的 FM 訊號，設定 \(f_m=1000Hz\)，\({\beta}=15\)，訊雜比 \(SNR=50dB\)，數位正交解調效果如下，解調後幅度放大係數 \(k=\overline{{\lvert}m(t){\rvert}}/\overline{{\lvert}\hat{m}(t){\rvert}}\approx2.99\)，使用這個係數放大解調訊號幅值，然後計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-k\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.0163\)。

數位正交解調的程式碼詳見 lpf_filter.m、demod_fm_method4.m 和 main_demodFM_example4.m。更改相干載波的初始相位為 \({\phi_0}=\pi/4,\pi/2\)，或者更改相干載波的中心頻率為 \(0.8f_c,1.2f_c\) 後，解調效果依然比較好，說明這種方法具有較好的抗載頻失配能力。

3. FM 模擬（MATLAB Communications Toolbox）

MATLAB 的 Communications Toolbox 中提供了 FM 調變函數 fmmod，高斯白噪聲函數 awgn，以及 FM 解調函數 fmdemod，可以很方便地完成 FM 訊號模擬。使用這三個函數實現上面 1.4 節中確知訊號 \(m(t)\) 的 FM 調變解調，調變後加噪聲的效果如下：

解調效果如下：

解調訊號與調變訊號波形基本重回，計算誤差，有：\(\sqrt{\sum{{\lvert}m(t_i)-\hat{m}(t_i){\rvert}^2}}/\sqrt{\sum{{\lvert}m(t_i){\rvert}^2}}\approx0.2397\)。程式碼詳見附錄 main_CommFM_example.m。

參考資料

[1] 樓才義,徐建良,楊小牛.軟體無線電原理與應用[M].電子工業出版社,2014.

[2] 樊昌信,曹麗娜.通訊原理.第7版[M].國防工業出版社,2012.

[3] 劉學勇.詳解MATLAB/Simulink通訊系統建模與模擬[M].電子工業出版社,2011.

[4] 王麗娜,王兵.衛星通訊系統.第2版[M].國防工業出版社,2014.

[5] 部落格園 - 兩種頻率調變(FM)方法的MATLAB實現。

[6] CSDN - 數位訊號處理基礎----FM的調變與解調。

附錄程式碼

附.1 檔案 mod_fm.m

function [ sig_fm, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A)
% MOD_FM        FM 調頻
% 輸入引數：
%       fc      載波中心頻率
%       beta    調頻指數/調變指數
%       fs      訊號取樣率
%       mt      調變訊號
%       t       取樣時間
%       W       基頻訊號頻寬（最高頻率）
%       A       載波恆定振幅
% 輸出引數：
%       sig_fm  調頻(FM)實訊號
%       deltaf  最大頻偏
% @author 木三百川

% 計算調頻靈敏度及最大頻偏
kf = beta*W/max(abs(mt));
deltaf = beta*W;

% 計算調變訊號積分
int_mt = cumtrapz(t,mt);

% 生成訊號
sig_fm = A*cos(2*pi*fc*t+2*pi*kf*int_mt); % FM調頻訊號

% 繪圖
nfft = length(sig_fm);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm));
subplot(3,2,1);
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('調變訊號m(t)');
subplot(3,2,2);
plot(freq, 10*log10(fftshift(abs(fft(mt,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('調變訊號m(t)雙邊幅度譜');

subplot(3,2,3);
plot(t(1:plot_length), int_mt(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('調變訊號m(t)積分');
subplot(3,2,4);
plot(freq, 10*log10(fftshift(abs(fft(int_mt,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('調變訊號m(t)積分雙邊幅度譜');

subplot(3,2,5);
plot(t(1:plot_length), sig_fm(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM調頻訊號s(t)');
subplot(3,2,6);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('FM調頻訊號s(t)雙邊幅度譜');

end

附.2 檔案 main_modFM_example1.m

clc;
clear;
close all;
% FM 調變模擬(調變訊號為確知訊號)
% @author 木三百川

% 調變引數
A = 1;                  % 載波恆定振幅
fm = 2500;              % 調變訊號引數
beta = 1e-6;            % 調頻指數/調變指數
fc = 20000;             % 載波頻率
fs = 8*fc;              % 取樣率
total_time = 2;         % 模擬時長，單位：秒

% 取樣時間
t = 0:1/fs:total_time-1/fs;

% 調變訊號為確知訊號
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;

% FM 調變
[ sig_fm, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大頻偏 deltaf = %.6f Hz.\n', deltaf);

附.3 檔案 main_modFM_example2.m

clc;
clear;
close all;
% FM 調變模擬(調變訊號為隨機訊號)
% @author 木三百川

% 調變引數
A = 1;                  % 載波恆定振幅
fH = 3000;          	% 基頻訊號頻寬
beta = 1e-6;            % 調頻指數/調變指數
fc = 20000;             % 載波頻率
fs = 8*fc;              % 取樣率
total_time = 2;         % 模擬時長，單位：秒

% 取樣時間
t = 0:1/fs:total_time-1/fs;

% 調變訊號為隨機訊號
mt = randn(size(t));
b = fir1(512, fH/(fs/2), 'low');
mt = filter(b,1,mt);
mt = mt - mean(mt);
W = fH;

% FM 調變
[ sig_fm, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大頻偏 deltaf = %.6f Hz.\n', deltaf);

附.4 檔案 lpf_filter.m

function sig_lpf = lpf_filter(sig_data, cutfre)
% LPF_FILTER    自定義理想低通濾波器
% 輸入引數：
%       sig_data        待濾波資料
%       cutfre          截止頻率，範圍 (0,1)
% 輸出引數：
%       sig_lpf         低通濾波結果
% @author 木三百川

nfft = length(sig_data);
lidx = round(nfft/2-cutfre*nfft/2);
ridx = nfft - lidx;
sig_fft_lpf = fftshift(fft(sig_data));
sig_fft_lpf([1:lidx,ridx:nfft]) = 0;
sig_lpf = real(ifft(fftshift(sig_fft_lpf)));

end

附.5 檔案 demod_fm_method1.m

function [ sig_fm_demod ] = demod_fm_method1(sig_fm_receive, fc, fs, t)
% DEMOD_FM_METHOD1        FM 非相干解調（鑑頻器）
% 輸入引數：
%       sig_fm_receive      FM 接收訊號，行向量
%       fc                  載波中心頻率
%       fs                  訊號取樣率
%       t                   取樣時間
% 輸出引數：
%       sig_fm_demod        解調結果，與 sig_fm_receive 等長
% @author 木三百川

% 第一步：求微分
sig_dfm = [diff(sig_fm_receive),0];

% 第二步：全波整流
sig_fm_abs = abs(sig_dfm);

% 第三步：低通濾波
sig_fm_lpf = lpf_filter(sig_fm_abs, fc/2/(fs/2));

% 第四步：去除直流分量
sig_fm_demod = sig_fm_lpf - mean(sig_fm_lpf);

% 繪圖
nfft = length(sig_fm_abs);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_abs));
subplot(3,2,1);
plot(t(1:plot_length), sig_fm_abs(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('全波整流結果');
subplot(3,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_abs,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('全波整流結果雙邊幅度譜');

subplot(3,2,3);
plot(t(1:plot_length), sig_fm_lpf(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('低通濾波結果');
subplot(3,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_lpf,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('低通濾波結果雙邊幅度譜');

subplot(3,2,5);
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('（去除直流）解調結果');
subplot(3,2,6);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_demod,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('（去除直流）解調結果雙邊幅度譜');

end

附.6 檔案 main_demodFM_example1.m

clc;
clear;
close all;
% FM 解調模擬(調變訊號為確知訊號，非相干解調)
% @author 木三百川

% 調變引數
A = 1;                  % 載波恆定振幅
fm = 1000;              % 調變訊號引數
beta = 15;              % 調頻指數/調變指數
fc = 20000;             % 載波頻率
fs = 8*fc;              % 取樣率
total_time = 2;         % 模擬時長，單位：秒

% 取樣時間
t = 0:1/fs:total_time-1/fs;

% 調變訊號為確知訊號
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;

% FM 調變
[ sig_fm_send, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大頻偏 deltaf = %.6f Hz.\n', deltaf);

% 加噪聲
snr = 50;               % 訊雜比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');

% 非相干解調
[ sig_fm_demod ] = demod_fm_method1(sig_fm_receive, fc, fs, t);

% 繪圖
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收訊號');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('FM接收訊號雙邊幅度譜');

coef = mean(abs(mt))/mean(abs(sig_fm_demod));
fprintf('norm(調變訊號 - %.2f * 解調訊號)/norm(調變訊號) = %.4f.\n', coef, norm(mt-coef*sig_fm_demod)/norm(mt));

figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), coef*sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解調效果');
legend('調變訊號','解調訊號(放大後)');

附.7 檔案 demod_fm_method2.m

function [ sig_fm_demod ] = demod_fm_method2(sig_fm_receive, fs, t)
% DEMOD_FM_METHOD2        FM 非相干解調（鑑頻器 - 希爾伯特變換法）
% 輸入引數：
%       sig_fm_receive      FM 接收訊號，行向量
%       fs                  訊號取樣率
%       t                   取樣時間
% 輸出引數：
%       sig_fm_demod        解調結果，與 sig_fm_receive 等長
% @author 木三百川

% 第一步：求微分
sig_dfm = [diff(sig_fm_receive),0];

% 第二步：計算訊號包絡
sig_fm_envelope = abs(hilbert(sig_dfm));

% 第三步：去除直流分量
sig_fm_demod = sig_fm_envelope - mean(sig_fm_envelope);

% 繪圖
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(3,2,1);
plot(t(1:plot_length), sig_dfm(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('求微分結果');
subplot(3,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_dfm,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('求微分結果雙邊幅度譜');

subplot(3,2,3);
plot(t(1:plot_length), sig_fm_envelope(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('計算訊號包絡結果');
subplot(3,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_envelope,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('計算訊號包絡結果雙邊幅度譜');

subplot(3,2,5);
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('（去除直流）解調結果');
subplot(3,2,6);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_demod,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('（去除直流）解調結果雙邊幅度譜');

end

附.8 檔案 main_demodFM_example2.m

clc;
clear;
close all;
% FM 解調模擬(調變訊號為確知訊號，非相干解調/希爾伯特變換法)
% @author 木三百川

% 調變引數
A = 1;                  % 載波恆定振幅
fm = 1000;              % 調變訊號引數
beta = 15;              % 調頻指數/調變指數
fc = 20000;             % 載波頻率
fs = 8*fc;              % 取樣率
total_time = 2;         % 模擬時長，單位：秒

% 取樣時間
t = 0:1/fs:total_time-1/fs;

% 調變訊號為確知訊號
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;

% FM 調變
[ sig_fm_send, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大頻偏 deltaf = %.6f Hz.\n', deltaf);

% 加噪聲
snr = 50;               % 訊雜比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');

% 非相干解調
[ sig_fm_demod ] = demod_fm_method2(sig_fm_receive, fs, t);

% 繪圖
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收訊號');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('FM接收訊號雙邊幅度譜');

coef = mean(abs(mt))/mean(abs(sig_fm_demod));
fprintf('norm(調變訊號 - %.2f * 解調訊號)/norm(調變訊號) = %.4f.\n', coef, norm(mt-coef*sig_fm_demod)/norm(mt));

figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), coef*sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解調效果');
legend('調變訊號','解調訊號(放大後)');

附.9 檔案 demod_fm_method3.m

function [ sig_fm_demod ] = demod_fm_method3(sig_fm_receive, fc, fs, t, phi0)
% DEMOD_FM_METHOD3        FM 相干解調
% 輸入引數：
%       sig_fm_receive      FM 接收訊號，行向量
%       fc                  載波中心頻率
%       fs                  訊號取樣率
%       t                   取樣時間
%       phi0                相干載波初始相位
% 輸出引數：
%       sig_fm_demod        解調結果，與 sig_fm_receive 等長
% @author 木三百川

% 第一步：乘以相干載波
sig_fm_ct = -sig_fm_receive.*sin(2*pi*fc*t+phi0);

% 第二步：低通濾波
sig_fm_lpf = lpf_filter(sig_fm_ct, fc/(fs/2));

% 第三步：求微分
sig_fm_demod = [diff(sig_fm_lpf),0]*fs;

% 繪圖
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(3,2,1);
plot(t(1:plot_length), sig_fm_ct(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('乘以相干載波結果');
subplot(3,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_ct,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('乘以相干載波結果雙邊幅度譜');

subplot(3,2,3);
plot(t(1:plot_length), sig_fm_lpf(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('低通濾波結果');
subplot(3,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_lpf,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('低通濾波結果雙邊幅度譜');

subplot(3,2,5);
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('（求微分）解調結果');
subplot(3,2,6);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_demod,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('（求微分）解調結果雙邊幅度譜');

end

附.10 檔案 main_demodFM_example3.m

clc;
clear;
close all;
% FM 解調模擬(調變訊號為確知訊號，相干解調)
% @author 木三百川

% 調變引數
A = 1;                  % 載波恆定振幅
fm = 2500;              % 調變訊號引數
beta = 1e-6;            % 調頻指數/調變指數
fc = 20000;             % 載波頻率
fs = 8*fc;              % 取樣率
total_time = 2;         % 模擬時長，單位：秒

% 取樣時間
t = 0:1/fs:total_time-1/fs;

% 調變訊號為確知訊號
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;

% FM 調變
[ sig_fm_send, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大頻偏 deltaf = %.6f Hz.\n', deltaf);

% 加噪聲
snr = 150;               % 訊雜比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');

% 相干解調
ini_phase = 0;
[ sig_fm_demod ] = demod_fm_method3(sig_fm_receive, fc, fs, t, ini_phase);

% 繪圖
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收訊號');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('FM接收訊號雙邊幅度譜');

coef = mean(abs(mt))/mean(abs(sig_fm_demod));
fprintf('norm(調變訊號 - %.2f * 解調訊號)/norm(調變訊號) = %.4f.\n', coef, norm(mt-coef*sig_fm_demod)/norm(mt));

figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), coef*sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解調效果');
legend('調變訊號','解調訊號(放大後)');

附.11 檔案 demod_fm_method4.m

function [ sig_fm_demod ] = demod_fm_method4(sig_fm_receive, fc, fs, t, phi0)
% DEMOD_FM_METHOD4        FM 數位正交解調/相干解調
% 輸入引數：
%       sig_fm_receive      FM 接收訊號，行向量
%       fc                  載波中心頻率
%       fs                  訊號取樣率
%       t                   取樣時間
%       phi0                相干載波初始相位
% 輸出引數：
%       sig_fm_demod        解調結果，與 sig_fm_receive 等長
% @author 木三百川

% 第一步：乘以正交相干載波
sig_fm_i = sig_fm_receive.*cos(2*pi*fc*t+phi0);
sig_fm_q = -sig_fm_receive.*sin(2*pi*fc*t+phi0);

% 第二步：低通濾波
sig_fm_i_lpf = lpf_filter(sig_fm_i, fc/(fs/2));
sig_fm_q_lpf = lpf_filter(sig_fm_q, fc/(fs/2));

% 第三步：計算相位
Pt = unwrap(atan2(sig_fm_q_lpf, sig_fm_i_lpf));

% 第四步：計算差分
sig_fm_demod = [diff(Pt),0];

% % 合併第三步與第四步：反正切近似
% sig_fm_demod = (sig_fm_i_lpf(1:end-1).*sig_fm_q_lpf(2:end)-sig_fm_i_lpf(2:end).* ...
%     sig_fm_q_lpf(1:end-1))./(sig_fm_i_lpf(2:end).^2+sig_fm_q_lpf(2:end).^2);
% sig_fm_demod = [sig_fm_demod, sig_fm_demod(end)];

% 繪圖
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(2,2,1);
plot(t(1:plot_length), sig_fm_i(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('乘以正交相干載波 I 路結果');
subplot(2,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_i,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('乘以正交相干載波 I 路結果雙邊幅度譜');
subplot(2,2,3);
plot(t(1:plot_length), sig_fm_q(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('乘以正交相干載波 Q 路結果');
subplot(2,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_q,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('乘以正交相干載波 Q 路結果雙邊幅度譜');

figure;set(gcf,'color','w');
subplot(2,2,1);
plot(t(1:plot_length), sig_fm_i_lpf(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('低通濾波 I 路結果');
subplot(2,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_i_lpf,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('低通濾波 I 路結果雙邊幅度譜');
subplot(2,2,3);
plot(t(1:plot_length), sig_fm_q_lpf(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('低通濾波 Q 路結果');
subplot(2,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_q_lpf,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('低通濾波 Q 路結果雙邊幅度譜');

figure;set(gcf,'color','w');
subplot(2,2,1);
plot(t(1:plot_length), Pt(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('計算相位結果');
subplot(2,2,2);
plot(freq, 10*log10(fftshift(abs(fft(Pt,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('計算相位結果雙邊幅度譜');
subplot(2,2,3);
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('（計算差分）解調結果');
subplot(2,2,4);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_demod,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('（計算差分）解調結果雙邊幅度譜');

end

附.12 檔案 main_demodFM_example4.m

clc;
clear;
close all;
% FM 解調模擬(調變訊號為確知訊號，數位正交解調)
% @author 木三百川

% 調變引數
A = 1;                  % 載波恆定振幅
fm = 1000;              % 調變訊號引數
beta = 15;              % 調頻指數/調變指數
fc = 20000;             % 載波頻率
fs = 8*fc;              % 取樣率
total_time = 2;         % 模擬時長，單位：秒

% 取樣時間
t = 0:1/fs:total_time-1/fs;

% 調變訊號為確知訊號
mt = sin(2*pi*fm*t)+cos(pi*fm*t);
W = fm;

% FM 調變
[ sig_fm_send, deltaf ] = mod_fm(fc, beta, fs, mt, t, W, A);
fprintf('最大頻偏 deltaf = %.6f Hz.\n', deltaf);

% 加噪聲
snr = 50;               % 訊雜比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');

% 相干解調
ini_phase = 0;
[ sig_fm_demod ] = demod_fm_method4(sig_fm_receive, fc, fs, t, ini_phase);

% 繪圖
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收訊號');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('FM接收訊號雙邊幅度譜');

coef = mean(abs(mt))/mean(abs(sig_fm_demod));
fprintf('norm(調變訊號 - %.2f * 解調訊號)/norm(調變訊號) = %.4f.\n', coef, norm(mt-coef*sig_fm_demod)/norm(mt));

figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), coef*sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解調效果');
legend('調變訊號','解調訊號(放大後)');

附.13 檔案 main_CommFM_example.m

clc;
clear;
close all;
% FM 調變解調模擬(使用Communications Toolbox工具箱)
% @author 木三百川

% 調變引數
A = 1;                  % 載波恆定振幅
fm = 2500;              % 調變訊號引數
beta = 1e-6;            % 調頻指數/調變指數
fc = 20000;             % 載波頻率
fs = 8*fc;              % 取樣率
total_time = 2;         % 模擬時長，單位：秒

% 取樣時間
t = 0:1/fs:total_time-1/fs;

% 調變訊號為確知訊號
mt = sin(2*pi*fm*t)+cos(pi*fm*t);

% FM 調變
freqdev = beta*fm;
ini_phase = 0;
sig_fm_send = A*fmmod(mt, fc, fs, freqdev, ini_phase);

% 加噪聲
snr = 150;               % 訊雜比
sig_fm_receive = awgn(sig_fm_send, snr, 'measured');

% FM 解調
[ sig_fm_demod ] = fmdemod(sig_fm_receive, fc, fs, freqdev, ini_phase);

% 繪圖
nfft = length(sig_fm_receive);
freq = (-nfft/2:nfft/2-1).'*(fs/nfft);
figure;set(gcf,'color','w');
plot_length = min(500, length(sig_fm_receive));
subplot(1,2,1);
plot(t(1:plot_length), sig_fm_receive(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('FM接收訊號');
subplot(1,2,2);
plot(freq, 10*log10(fftshift(abs(fft(sig_fm_receive,nfft)/nfft))+eps));xlim([freq(1),freq(end)]);
xlabel('頻率/hz');ylabel('幅度/dB');title('FM接收訊號雙邊幅度譜');

figure;set(gcf,'color','w');
plot(t(1:plot_length), mt(1:plot_length));xlim([t(1),t(plot_length)]);
hold on;
plot(t(1:plot_length), sig_fm_demod(1:plot_length));xlim([t(1),t(plot_length)]);
xlabel('t/s');ylabel('幅度');title('解調效果');
legend('調變訊號','解調訊號');

fprintf('norm(調變訊號 - 解調訊號)/norm(調變訊號) = %.4f.\n', norm(mt-sig_fm_demod)/norm(mt));