線性調頻信號（chirp信號）

顧名思義，該信號的頻率隨著時間線性變換，其復數(shù)表達形式如下：
$$s(t)=e^{2j\pi (f_{0}t+0.5\mu t^2)}$$
根據(jù)歐拉公式，其相位項為$2\pi (f_{0}t+0.5\mu t^2))$。信號的角頻率是相位對時間的導數(shù)。對相位求導后有$2\pi (f_0+\mu t)$，顯然，$f_0$表示線性調頻起始頻率，$\mu$表示調頻斜率。

下面設一個線性調頻信號，起始頻率$f_0=24GHz$,調頻帶寬$B=10GHz$，采樣頻率為$f_s=240GHz$，采樣時寬$T=100ns$，據(jù)此可以計算調頻斜率為$\mu =B/T=10^{17}Hz/s$（注：參數(shù)僅供參考，并未考慮實際可行性）

對該信號添加一個信噪比為$5$dB的高斯白噪聲，可以得到其時域圖如下所示：

上方的圖是前1e-9秒的信號，下放的圖是最后1e-9秒的信號，可以看出信號變得“緊密”了一些，這是因為頻率從$24GHz$線性上升到了$34GHz$。

利用fft觀察雙邊頻譜如下圖

雖然噪聲很大，但是也可以大致看出，在$24GHz$到$34GHz$的頻帶上幅度較高。
這張圖還可以說明，一般的傅里葉變換對chirp信號的檢測能力不足，噪聲稍大一些就會淹沒信號。下面利用分數(shù)階傅里葉變換處理chirp信號。

分數(shù)階傅里葉變換檢測chirp信號

線性調頻信號作為一種典型的非平穩(wěn)信號，具有大時寬帶寬積的特殊優(yōu)勢，廣泛應用于雷達、通信、地質探測和聲吶信號處理等研究領域。因此，研究線性調頻信號具有十分重要的意義。分數(shù)階傅里葉變換實質上是一種線性變換,它不僅可以理解為chirp基分解，還沒有交叉干擾的問題。所以，分數(shù)階傅里葉變換特別適合用來處理chirp類信號。[1]

我的理解：就像平面任意一個矢量可以分解為x,y兩個基向量一樣，任意一個線性調頻信號也可以分解成兩個chirp基。只要選擇對了合適的階數(shù)，分數(shù)階傅里葉變換就可以在把這個chirp信號的能量集中在這個方向上的基向量上，形成一個沖激脈沖（易于檢測），而一般的傅里葉變換通常并不在這個“最優(yōu)方向”上，所以能量散開，具有了一定的帶寬，不利于檢測。
（以上想法僅供參考，不保證正確性）

分數(shù)階傅里葉變換的原理及其快速算法（Ozakats算法）的MATLAB實現(xiàn)在我的另一篇博客中：分數(shù)階傅里葉變換（FrFT）詳細原理與matlab代碼實現(xiàn)

下面從最基礎的分數(shù)階傅里葉變換公式入手[2]：

式中的$B_a(x,x^{\prime })$表示卷積核，$f(x)$表示信號，$x^{\prime }$表示與$x$不同的變量而不是導數(shù)。
將上述線性調頻信號的公式代入，合并指數(shù)項，可以得到積分號里面的表達式為：
$$\begin{gathered} B_a(x,x^{\prime })f(x^{\prime })=A_{?}\exp (i\pi (x^2\cot ?)?2x{x}^{\prime }\csc ?+x^{\prime 2}\cot ?)?\exp (i\pi \mu x^{\prime 2}+2i\pi f_0{x}^{\prime }) \\ =A_{?}\exp (i\pi (x^2\cot ?)+2x^{\prime }(f_0?x\csc ?)+x^{\prime 2}(\mu +\cot ?)) \end{gathered}$$
該積分是對$x^{\prime }$變量積分，于是著重觀察含$x^{\prime }$的項，包括一次項$2x^{\prime }(f_0?x\csc ?)$和二次項$x^{\prime 2}(\mu +\cot ?)$.

調整階數(shù)$a$，使得$\mu =?\cot ?$，消去二次項。

在沒有二次項的條件下，再假設$x=f_0/\csc ?$，消去一次項，于是積分內部沒有$x^{\prime }$，可以去掉積分號，變換結果為$A_{?}\exp (i\pi x^2\cot ?)$，代入$x=f_0/\csc ?$，顯然這是一個非零項。
如果$x\ne f_0/\csc ?$，就會出現(xiàn)一個一次項的積分。由于$\sin t$和$\cos t$的周期性，其無窮積分等于零，再根據(jù)歐拉公式，有下式成立:
$${\int }_{?\infty }^{\infty }{e}^{it}dt=0$$
因此一次項的存在會使整個積分為0，換句話說，該積分只有在$x=f_0/\csc ?$處有值，而在其它地方為0，類似于沖激函數(shù)。
在二次項存在的情況下，該積分就沒有上述性質，能量散開，不易檢測。

參考文獻[1]總結了利用分數(shù)階傅里葉變換檢測chirp信號的步驟以及參數(shù)估計的公式：

其中${\alpha }_0$在上文中指$?$。整個問題變成了一個二維最優(yōu)化問題，優(yōu)化目標是最大化分數(shù)階傅里葉變換函數(shù)的值，參數(shù)為階數(shù)$a$和頻率$x$。

注意，在實際操作中發(fā)現(xiàn)上述參數(shù)估計方程并不正確。進一步查閱文獻并參考了他人代碼后發(fā)現(xiàn)，上述參數(shù)估計方程是適用于連續(xù)函數(shù)的，在實際的計算機中都是離散信號，因此${\mu }_0$需要使用歸一化調頻斜率$k=\frac{B}{f_s}$，且$x\csc ?$實際是中心頻率而不是起始頻率。

---

2022.10.3學習更新：
此處得出的確實是中心頻率，原因是本文的信號是[0,T]范圍內的，而上述推導過程包括參數(shù)歸一化是在[-T/2,T/2]范圍內的，當信號時間處于這個范圍內，估計出來的頻率確實是初始頻率。
有前輩也有這個疑惑，參見matlab論壇https://www.ilovematlab.cn/thread-277599-1-1.html

---

仿真

在例如汽車防撞雷達等一些場景中，我們提前知道要檢測的chirp信號參數(shù)，可以直接根據(jù)參數(shù)估計分數(shù)階傅里葉分解的階數(shù)，進而直接檢測目標信號。
假設我們使用上文中的chirp參數(shù)，首先求解歸一化調頻斜率$k=\frac{B}{f_s}=0.0417$
代入參數(shù)估計公式中，計算旋轉角度$?=arccot(?k)=?1.5292$，然后計算階次$a=?/\frac{\pi }{2}=?0.9735$
負數(shù)不影響結果，實際上分數(shù)階傅里葉變換周期為4，而對于chirp信號檢測來說，$\cot (?)=\cot (?+\pi )，$周期進一步縮小為2。

對chirp信號做階數(shù)為-0.9735的分數(shù)階傅里葉變換，得到如下時頻域圖：

由圖可以看出，信號在2.8976e10的位置出現(xiàn)了一個非常大的峰值，相比于fft的頻域圖像，顯然這個峰值更容易檢測到。因此使用分數(shù)階傅里葉變換比一般的傅里葉變換更抗噪聲干擾。

令$x=2.8976e10$，估計信號的中心頻率為$x\csc ?=29Ghz$，因此可以推算出起始頻率為24Ghz，與預設相符。

代碼

close all fs=24e10;%采樣頻率 T=1e-7;%時寬 B=10e9;%帶寬 mu=B/T;%調頻率 n=round(T*fs);%采樣點個數(shù) t=linspace(0,T,n); f0=24e9;%起始頻率 s=exp(2j*pi*(f0*t+0.5*mu*t.^2)); SNR=5;%信噪比 s=awgn(s,SNR);%添加一定信噪比的高斯白噪聲figure subplot(211) plot(t,real(s))%時域圖 title("調頻信號時域") xlabel("t/s") xlim([0,1e-9]) grid onsubplot(212) plot(t,real(s))%時域圖 title("調頻信號時域") xlabel("t/s") xlim([T-1e-9,T]) grid onfigure S=fftshift(fft(real(s))./n); f=linspace(-fs/2,fs/2-1,n);%頻域橫坐標，注意奈奎斯特采樣定理，最大原信號最大頻率不超過采樣頻率的一半 plot(f,abs(S))%頻域圖 title("發(fā)射信號頻域") xlim([-50e9,50e9]) grid onk=B/fs;%歸一化調頻斜率 a=acot(-k)/(pi/2);%以歸一化調頻斜率估計階次 % a=pi/2+atan(n-1)*mu/fs^2; % for a=0.9:0.01:1 figure S=myfrft(real(s),a); f=linspace(-fs/2,fs/2-1,n);%頻域橫坐標，注意奈奎斯特采樣定理，最大原信號最大頻率不超過采樣頻率的一半 plot(f,abs(S))%頻域圖 title("發(fā)射信號時頻域") % xlim([0,fs/2-1]) title("a="+num2str(a)) grid on % endfh_=abs(f(abs(S)==max(abs(S)))*csc(a*pi/2));%估計的中心頻率 f0_=fh-B/2;%估計的起始頻率

參考文獻

[1] 渠瑩,楊俊. 基于分數(shù)階傅里葉變換的LFM信號參數(shù)估計[J]. 物聯(lián)網技術,2017,7(11):30-32. DOI:10.16667/j.issn.2095-1302.2017.11.006.
[2] OZAKTAS H.M., ARIKAN O… Digital computation of the fractional Fourier transform[J]. IEEE Transactions on Signal Processing: A publication of the IEEE Signal Processing Society,1996,44(9):2141-2150.

---

注：本人系初學者，以上僅為個人學習筆記，可能存在錯誤，有任何問題歡迎在評論區(qū)討論。本文將會隨著學習理解的加深而繼續(xù)修改。

總結

以上是生活随笔為你收集整理的基于分数阶傅里叶变换的chirp信号检测与参数估计（原理附代码）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。