為什么要寫這一篇文章呢？其實最開始我只是想學習下鎖相環，然后用MATLAB做一下仿真什么的。但我認為最好的學習方法是找一個簡單的應用，這應用里要用到鎖相環的相關知識。我想到了頻率調制在同步解調時可能會用到它，於是去查找了相關資料，實現了部分內容並寫下了這篇文章，文章中為了完整的展現整個過程，索性將完整的FM調制解調都總結了一下。為了說明代碼的計算原理，也進行了一些相關知識的補充。
一: 頻率調制
1.1 要實現什么
在通信原理這門課中，我們學習知道了頻率調制FM的目的，是用載波的頻率大小來表示基帶信號幅度的大小。比如下面這幅圖，可以看到我們想要實現得到一個已調波，基帶信號幅度上升時，已調波頻率增加；基帶信號幅度下降時，已調波頻率減小。
1.2 原理與實現
那么如何用基帶信號的幅度來控制已調波的頻率呢？那我們就人為的把他倆扯上關系吧。下面是一個把他倆扯上關系的實現過程：
首先我們假設基帶信號為\(m(t)\)，載波頻率為\(f_c\)，載波振幅為\(A_c\)。我們倒着推，直接假設已調波\(S_{FM}(t)\)為：
瞬時頻率為\(f_k = 2\pi f_c + \frac{d\phi(t)}{dt}\)
瞬時頻偏為\(\frac{d\phi(t)}{dt}\)
把瞬時頻偏和信號的幅度扯上關系，便滿足了這一節剛開始我們的想法。設下面關系式：
其中\(K_f\)為頻偏常數，表示調頻器的調頻靈敏度。
將(2)代回(1)里面，得到已調波的真正表達式為：
更進一步，假設基帶信號為余弦波\(m(t) = cos(2 \pi f_m t)\)，\(A_c = 1\)，已調波的表達式為：
順利推導出了已調波的表達式，調制部分差不多就結束了，我們直接根據表達式寫代碼就好了。
Fs = 2000; % 采樣頻率
Ts = 1 / Fs;
L = 2000; % 取2000采樣點，1s
t = (0 : L - 1) * Ts;
kf = 45; % 調頻靈敏度
fm = 4; % 基帶信號頻率
fc = 50; % 載波頻率
mt = cos(2 * pi * fm * t); % 基帶信號
st = cos(2 * pi * fc * t + kf * sin(2 * pi * fm * t) / fm); % 由公式(5)寫出已調FM信號
% st = cos(2 * pi * fc * t + 2 * pi * kf * cumsum(mt) * Ts); % 也可由公式(4)寫出已調FM信號
下圖是基帶信號和已調信號波形圖：
二: 非相干解調
現在基帶信號已經調制完畢了。我們不考慮噪聲的問題，之后該如何從已調信號恢復我們想要的基帶信號呢？首先我們使用非相干解調的方法，顧名思義，使用非相干解調的話，我們就不用進行把已調信號和本地載波相乘的操作，也不用使用鎖相環。但我們需要對已調信號做微分和使用一點點希爾伯特變換相關的知識，下面將進行詳細講解。
2.1 要實現什么
我們要想辦法從已調信號中恢復出\(m(t)\)來，對公式\((3)\)的\(S_{FM}(t)\)求微分，得到：
一般假設\(A_c = 1\)，現在我們就要思考了，要是我們能將\(S_d(t)\)的包絡\(2 \pi f_c + K_f m(t)\)通過某種方式提取出來，那么我們不就能輕易的恢復\(m(t)\)了嗎？沒錯，接下來我們就要用希爾伯特變換和解析信號的性質來提取包絡了！
等下！我們先來看看結論，如圖可知如果我們能夠提取出包絡，便能很容易的恢復基帶信號了，你看包絡和基帶信號長得多像，只是振幅小了些，增加了直流分量而已。
2.2 希爾伯特變換(包絡檢波原理)
這里補充一些希爾伯特變換的相關知識。首先為什么研究學者要引入希爾伯特變換這個玩意兒呢？它其實是有實際的現實意義的。
一般我們在計算某一信號的頻譜時，我們會發現信號除了正的頻率外，還有對稱的另一半負頻率，比如上面公式\((6)\)的微分信號\(S_d(t)\)的頻譜圖如下所示：
我們設某一信號的頻譜為\(X(f)\)，如果我們想去掉負頻率部分，並且能量不變的話，我們可以將原頻譜的2倍乘以階躍函數\(U(f)\)，於是有：
上式是頻域的變化，那時域是如何變的呢，由傅里葉逆變換得到：
於是我們知道了，想要除去信號的負頻率部分，只需要用該信號與\(\frac{1}{\pi t}\)的卷積乘以i，再加上該信號即可。一般我們把\(\widetilde x(t)\)叫做\(x(t)\)的解析信號；把\(x(t) * \frac{1}{\pi t}\)叫做\(x(t)\)的希爾伯特變換。
上面就是希爾伯特變換的來源和作用了（這是附加內容，了解下就可以了），定義希爾伯特變換:
定義解析信號：
那么如何用解析信號來求包絡呢？
假設包絡為\(a(t)\)，載波為\(cos(2\pi f_0 t + \theta(t))\)，已調信號\(x(t) = a(t)cos(2\pi f_0 t + \theta(t))\)。
解析信號為：
取解析信號的模：
可以看到我們成功提取出了包絡來。其實整個過程就是我們構造解析信號，然后利用的希爾伯特變換的正交性，從已調信號中提取出包絡。
2.3 包絡檢波實現
由上面的原理我們已經知道了如何從已調信號中提取出包絡來，直接上實現代碼吧。代碼超簡單，就兩行QAQ。注意MATLAB中的hilbert函數是求解析信號而不是求希爾伯特變換，要得到希爾伯特變換的話由（10）可知，取其中的虛部就可以了。
sdt = [0 diff(st)]; % 先對已調信號做微分
rt = abs(hilbert(sdt)); % 再構造解析信號求包絡,得到非相干解調信號
從上圖可以看到我們已經基本恢復了基帶信號，剩下的只需做一下隔直流和把信號放大就可以了。注意：解調信號振幅為什么變小呢？因為在MATLAB做微分時，用函數diff()並不能對時間間隔取到無窮小。
三: 相干解調(同步解調)
在這一部分中，要考慮用另外一種FM解調方式了，就是利用鎖相環(PLL)進行相干解調。鎖相環的作用之一便是在信號解調方面的應用。對於鎖相環的原理和實現，參考了很多相關資料，在下面是我的一個粗淺的總結。
3.1 鎖相環功能與結構
鎖相環是一種相位反饋跟蹤系統。鑒相器輸出的相位誤差信號經過環路濾波器濾波后，作為壓控振盪器的控制信號，而壓控振盪器的輸出又反饋到鑒相器，在鑒相器中與輸入信號進行相位比較。壓控振盪器的輸出信號相位將跟蹤輸入信號的相位變化，這時壓控振盪器輸出信號的頻率與輸入信號頻率相等，而相位保持一個微小誤差。
PLL一般用於：
- 載波同步
- FM解調
這里我們直接按照FM解調的數學模型討論。框圖如下：
PLL主要包括3部分：
-
壓控振盪器VCO(Voltage Controlled Oscillator)：根據輸入信號大小的不同，輸出頻率不同的的正弦波。
-
乘法鑒相器PD(Phase Detector)：將已調信號\(s(t)\)與VCO的輸出信號\(r(t)\)相乘。
-
環路濾波器LF(Loop Filter)：低通濾波器，將PD輸出的高頻分量濾除。
3.2 要實現什么
鎖相環所要實現的是：給我們一個輸入信號\(s(t)\)，\(s(t)\)的瞬時相位可能是在不斷變化着的，但我們通過這樣一個鎖相環結構，使得VCO的輸出\(r(t)\)的瞬時相位跟隨和追蹤\(s(t)\)，當跟蹤成功時，這時\(r(t)\)基本上和\(s(t)\)一個樣，此時我們也能從環路中恢復出基帶信號\(m(t)\)了，具體如何恢復的見下面推導。
3.3 鎖相環各部分公式推導
3.3.1 VCO部分
我們由公式（3）已知輸入的FM已調波為(正弦余弦都可以，為了好說明，這里用正弦):
其中
一般假定VCO的固有振盪頻率為載波頻率\(f_c\)，設由本地VCO跟蹤產生的FM信號為：
其中
\(A_v\)為VCO載波幅度，\(K_v\)為VCO的頻率靈敏度。從公式的角度，鎖相環的目的就是不斷調整\(v(t)\)，使得\(\phi_2(t)\)盡可能等於\(\phi_1(t)\)
3.3.2 PD部分
將輸入的已調波\(s(t)\)與本地FM波相乘，會得到兩個分量。
3.3.3 LF部分
通過低通濾波器，濾去\(e(t)\)的高頻分量，輸出為\(v(t)\)。
\(K_m\)為鑒相器增益；\(\phi_e(t) = \phi_1(t) - \phi_2(t)\)，為相位誤差。
當\(\phi_e(t) = 0\)，就說PLL處於鎖定狀態；當\(\phi_e(t) < 1 rad\)，就說PLL處於接近鎖定狀態，並且可以近似運算\(sin(\phi_e(t)) \approx \phi_e(t)\)。所以(17)進一步寫為：
由於\(\phi_2(t)\)在不斷跟蹤\(\phi_1(t)\)，所以處於穩定狀態時，\(\phi_e(t) =0\)，\(\phi_2(t) = \phi_1(t)\)，由公式(13)(15)可得：
就能用\(v(t)\)恢復\(m(t)\)了。
3.3.4 低通濾波器原理與實現補充
在PLL中會用到低通濾波器，為了之后能看懂低通濾波的代碼，這里介紹一下IIR濾波器原理，如何使用濾波系數進行濾波，以及如何使用MATLAB的FDA工具獲取濾波系數，這里想多總結一點，不想看的話可以有個印象即可。
通常系統(線性時不變離散系統)輸入輸出之間的運算關系可由N階線性常系數差分方程來描述：
比如當\(M = 0\)且\(N = 1\)，系統為\(y[n] = b_0x[n] -a_1y[n - 1]\)。
同樣的，系統也可以由輸入\(x[n]\)與單位脈沖響應\(h[n]\)的卷積來描述：
兩邊取\(Z\)變換
其中\(H(z) = \frac{Y(z)}{X(z)}\)定義為系統函數，濾波器就是一種系統，也是由系統函數來刻畫的。
- 濾波器當前的輸出依賴於輸入和過去的輸出，如\((19)\)，稱為無限長單位脈沖響應濾波器（IIR）。
- 濾波器當前的輸出僅依賴於輸入，而不依賴過去的輸出，比如\((19)\)中\(a_i = 0\)的情況，稱為有限長單位脈沖響應濾波器（FIR）。
於是IIR數字濾波器可用\((19)\)的差分方程描述。其中\(b_i\)和\(a_i\)決定了每個輸入和輸出的貢獻大小，稱為濾波器系數，\(N\)為濾波器階數。
IIR數字濾波器也可以用系統函數描述，對\((19)\)兩邊做\(Z\)變換：
推出IIR的系統函數為：
我們一般用計算機輔助設計濾波器時，我們設定好相關參數，結果返回的就是\(b_i\)和\(a_i\)，然后我們再把輸入輸出經過\((19)\)過程，就完成了濾波。
但如果想知道濾波器系數\(b_i\)和\(a_i\)的具體計算過程的話，建議拿出手邊的 《數字信號處理教程》 ，翻到講解IIR濾波器設計的章節，過程涉及拉氏變換、\(Z\)變換、模擬濾波器映射到數字濾波器雙線性變換過程。在百度上是很難搜到計算步驟的，基本沒人寫這個。如果沒書的話可以參考巴特沃斯低通濾波器設計。
如下是用MATLAB輔助設計的步驟，就很容易了：
假設我們需要一個一階巴特沃斯低通濾波器，N = 1，采樣頻率Fs=40000Hz，截止頻率Fcut = 1000Hz。
令M = 1，根據公式\((19)\)，濾波器的差分方程為：
只要計算出濾波器系數\(b_0\)、\(b_1\)、\(a_1\)，我們就可以直接用差分方程計算濾波器輸出結果了，下面是如何使用MATLAB計算它們。
首先在MATLAB中打開Filter Design & Analysis，設置參數如下圖，然后在上圖標中選擇Filter Coefficients，再在Edit中選擇Convert to Single Section，然后看到的Numerator對應濾波器系數\(b_i\)，Denominator對應\(a_i\)。濾波系數都有了，接下來該做什么？運行一遍\((22)\)就完事兒了。
3.4 總結
這是最后的部分輸出波形圖：
從上圖可以看到，VCO對已調信號的跟蹤效果還是不錯的。
從上圖可以看到，鎖相環中LF的輸出v(t)確實如公式(18)所示，成功恢復了基帶信號。
四：完整代碼
這里與前面的部分代碼相比更改了下采樣頻率和載波頻率以及基帶信號頻率，還有已調波改成了sin，其余一樣。
clc, clear all, close all;
Fs = 40000; % 采樣頻率
Ts = 1 / Fs;
L = 5000;
t = (0 : L - 1) * Ts;
%% 頻率調制(FM)
kf = 2000; % 調頻靈敏度
fm = 500; % 基帶信號頻率
fc = 10000; % 載波頻率
mt = cos(2 * pi * fm * t); % 基帶信號
st = sin(2 * pi * fc * t + kf * sin(2 * pi * fm * t) / fm); % 由公式(5)寫出已調FM信號
% st = sin(2 * pi * fc * t + 2 * pi * kf * cumsum(mt) * Ts); % 也可由公式(4)寫出已調FM信號
figure;
subplot(2,1,1)
plot(t,mt); % 觀察基帶信號時域
title('基帶信號時域圖'); xlabel('時間(s)'); ylabel('幅度');
subplot(2,1,2)
plot(t,st); % 觀察已調信號時域
title('已調信號時域圖'); xlabel('時間(s)'); ylabel('幅度');
[ST,f1] = myfft_plot(st,Fs); % 觀察已調信號st的頻譜
figure;
plot(f1,abs(ST));
title('已調信號頻譜圖'); xlabel('頻率(Hz)'); ylabel('幅度');
%% 非相干解調(包絡檢波法)
sdt = [0 diff(st)]; % 先對已調信號做微分,公式(6)
rt = abs(hilbert(sdt)); % 再構造解析信號求包絡,得到非相干解調信號，公式(11)
% test = 2 * pi * fc + kf .* mt .\ 2000;
figure;
subplot(2,1,1);
plot(t,sdt,t,rt,'r');
legend('微分信號','信號包絡');
title('已調信號微分后時域圖'); xlabel('時間(s)'); ylabel('幅度');
subplot(2,1,2)
plot(t,mt,'r'); % 觀察已調信號時域
title('基帶信號時域圖'); xlabel('時間(s)'); ylabel('幅度');
figure;
[RT3,f4] = myfft_plot(sdt,Fs); % 觀察微分信號頻譜
subplot(2,1,1);
plot(f4,abs(RT3));
title('微分信號Sd(t)頻譜圖'); xlabel('頻率(Hz)'); ylabel('幅度');
[RT4,f5] = myfft_plot(hilbert(sdt),Fs); % 觀察微分信號單邊頻譜
subplot(2,1,2);
plot(f5,abs(RT4));
title('微分信號Sd(t)單邊頻譜圖'); xlabel('頻率(Hz)'); ylabel('幅度');
[RT,f2] = myfft_plot(rt,Fs); % 觀察解調信號頻譜
figure;
plot(f2,abs(RT));
title('解調信號頻譜圖'); xlabel('頻率(Hz)'); ylabel('幅度');
figure;
subplot(2,1,1);
plot(t,rt); % 觀察解調信號與基帶信號
title('解調信號'); xlabel('時間(s)'); ylabel('幅度');
subplot(2,1,2);
plot(t,mt);
title('基帶信號'); xlabel('時間(s)'); ylabel('幅度');
%% 相干解調法(同步解調)
% 用鎖相環同步
vco_phase = zeros(1,L); % 初始化vco相位
rt = zeros(1,L); % 初始化壓控振盪器vco輸出
et = zeros(1,L); % 初始化乘法鑒相器pd輸出
vt = zeros(1,L); % 初始化環路濾波器lf輸出
Av = 1; % vco輸出幅度
kv = 40000; % vco頻率靈敏度（）
km = 1; % pd增益
k0 = 1; % lf增益
rt(1) = Av * cos(vco_phase(1)); % vco第一次輸出
et(1) = km * rt(1) * st(1); % pd第一次輸出
b0 = 0.07295965726826667; % Fs = 40000，fcut = 1000的1階巴特沃斯低通濾波器系數，由FDA生成
b1 = 0.07295965726826667;
a1 = -0.8540806854634666;
vt(1) = k0 * (b0 * et(1)); % 由數字濾波器的差分方程計算lf第一次輸出，式(22)
for n = 2 : L; % 不斷計算vco的下一次輸出相位，最終會跟蹤已調信號的相位
vco_phase_change = 2 * pi * fc * Ts + kv * vt(n - 1) * Ts;
vco_phase(n) = vco_phase(n - 1) + vco_phase_change;
rt(n) = Av * cos(vco_phase(n)); % vco輸出（會跟蹤st的相位）
et(n) = km * rt(n) * st(n); % 乘法鑒相器輸出，式(16)
vt(n) = k0 * (b0 * et(n) + b1 * et(n - 1) - a1 * vt(n - 1));% 由數字濾波器的差分方程計算lf輸出，式(22)
end
figure;
subplot(2,1,1);
plot(t,st);
title('已調信號m(t)'); xlabel('時間(s)'); ylabel('幅度'); xlim([0 0.01]);
subplot(2,1,2);
plot(t,rt); % 觀察vco輸出信號與已調信號
title('vco輸出信號r(t)'); xlabel('時間(s)'); ylabel('幅度'); xlim([0 0.01]);
figure
plot(t,vt,t,mt)
legend('lf輸出信號v(t)','基帶信號m(t)');
xlim([0 0.05]);
添加一段myfft_plot函數的代碼如下：
%----------------------------------------------------
%介紹:計算信號的雙邊幅度頻譜和其橫坐標、並調整使得橫坐標中心頻率為0Hz
%
%輸入：x 輸入信號
% Fs 采樣頻率
%
%輸出：X 輸入信號的幅度頻譜（經過fftshift的）
% freq 輸入信號幅度頻譜的橫坐標
%----------------------------------------------------
function [X,freq]=myfft_plot(x,Fs)
N=length(x);
% 這部分計算頻譜的橫坐標，使得中心頻率為0Hz
if mod(N,2)==0
k=-N/2:N/2-1; % N為偶數
else
k=-(N-1)/2:(N-1)/2; % N為奇數
end
T=N/Fs;
freq=k/T; % 頻譜的橫坐標
% takes the fft of the signal, and adjusts the amplitude accordingly
X=fft(x)/N; % fft並歸一化
X=fftshift(X); % shifts the fft data so that it is centered
五: 參考文獻
[1] 吳茂. MATLAB R2016a 通信系統建模與仿真 28 個案例分析. 清華大學出版社, 2018.