基于分数阶傅里叶变换的chirp信号检测与参数估计（原理附代码）

线性调频信号（chirp信号）

顾名思义，该信号的频率随着时间线性变换，其复数表达形式如下：

s ( t ) = e 2 j π ( f 0 t + 0.5 μ t 2 ) s(t)=e^{2j\pi(f_0 t+ 0.5\mu t^2)} s(t)=e2jπ(f0​t+0.5μt2)

根据欧拉公式，其相位项为 2 π ( f 0 t + 0.5 μ t 2 ) ) 2\pi(f_0 t+ 0.5\mu t^2)) 2π(f0​t+0.5μt2))。信号的角频率是相位对时间的导数。对相位求导后有 2 π ( f 0 + μ t ) 2\pi(f_0+\mu t) 2π(f0​+μt)，显然， f 0 f_0 f0​表示线性调频起始频率， μ \mu μ表示调频斜率。

下面设一个线性调频信号，起始频率 f 0 = 24 G H z f_0=24GHz f0​=24GHz,调频带宽 B = 10 G H z B=10GHz B=10GHz，采样频率为 f s = 240 G H z f_s=240GHz fs​=240GHz，采样时宽 T = 100 n s T=100ns T=100ns，据此可以计算调频斜率为 μ = B / T = 1 0 17 H z / s \mu=B/T=10^{17} Hz/s μ=B/T=1017Hz/s（注：参数仅供参考，并未考虑实际可行性）

对该信号添加一个信噪比为 5 5 5dB的高斯白噪声，可以得到其时域图如下所示：

上方的图是前1e-9秒的信号，下放的图是最后1e-9秒的信号，可以看出信号变得“紧密”了一些，这是因为频率从 24 G H z 24GHz 24GHz线性上升到了 34 G H z 34GHz 34GHz。

利用fft观察双边频谱如下图

虽然噪声很大，但是也可以大致看出，在 24 G H z 24GHz 24GHz到 34 G H z 34GHz 34GHz的频带上幅度较高。

这张图还可以说明，一般的傅里叶变换对chirp信号的检测能力不足，噪声稍大一些就会淹没信号。下面利用分数阶傅里叶变换处理chirp信号。

分数阶傅里叶变换检测chirp信号

线性调频信号作为一种典型的非平稳信号，具有大时宽带宽积的特殊优势，广泛应用于雷达、通信、地质探测和声呐信号处理等研究领域。因此，研究线性调频信号具有十分重要的意义。分数阶傅里叶变换实质上是一种线性变换,它不仅可以理解为chirp基分解，还没有交叉干扰的问题。所以，分数阶傅里叶变换特别适合用来处理chirp类信号。[1]

我的理解：就像平面任意一个矢量可以分解为x,y两个基向量一样，任意一个线性调频信号也可以分解成两个chirp基。只要选择对了合适的阶数，分数阶傅里叶变换就可以在把这个chirp信号的能量集中在这个方向上的基向量上，形成一个冲激脉冲（易于检测），而一般的傅里叶变换通常并不在这个“最优方向”上，所以能量散开，具有了一定的带宽，不利于检测。

（以上想法仅供参考，不保证正确性）

分数阶傅里叶变换的原理及其快速算法（Ozakats算法）的MATLAB实现在我的另一篇博客中：分数阶傅里叶变换（FrFT）详细原理与matlab代码实现

下面从最基础的分数阶傅里叶变换公式入手[2]：

式中的 B a ( x , x ′ ) B_a(x,x') Ba​(x,x′)表示卷积核， f ( x ) f(x) f(x)表示信号， x ′ x' x′表示与 x x x不同的变量而不是导数。

将上述线性调频信号的公式代入，合并指数项，可以得到积分号里面的表达式为：

B a ( x , x ′ ) f ( x ′ ) = A ϕ exp ⁡ ( i π ( x 2 cot ⁡ ϕ ) − 2 x x ′ csc ⁡ ϕ + x ′ 2 cot ⁡ ϕ ) ⋅ exp ⁡ ( i π μ x ′ 2 + 2 i π f 0 x ′ ) = A ϕ exp ⁡ ( i π ( x 2 cot ⁡ ϕ ) + 2 x ′ ( f 0 − x csc ⁡ ϕ ) + x ′ 2 ( μ + cot ⁡ ϕ ) ) B_a(x,x')f(x')=A_\phi \exp(i\pi(x^2\cot\phi)-2xx'\csc\phi+x'^2\cot\phi)\cdot\exp(i\pi\mu x'^2+2i\pi f_0x')\\=A_\phi \exp(i\pi(x^2\cot\phi)+2x'(f_0-x\csc\phi)+x'^2(\mu+\cot\phi)) Ba​(x,x′)f(x′)=Aϕ​exp(iπ(x2cotϕ)−2xx′cscϕ+x′2cotϕ)⋅exp(iπμx′2+2iπf0​x′)=Aϕ​exp(iπ(x2cotϕ)+2x′(f0​−xcscϕ)+x′2(μ+cotϕ))

该积分是对 x ′ x' x′变量积分，于是着重观察含 x ′ x' x′的项，包括一次项 2 x ′ ( f 0 − x csc ⁡ ϕ ) 2x'(f_0-x\csc\phi) 2x′(f0​−xcscϕ)和二次项 x ′ 2 ( μ + cot ⁡ ϕ ) x'^2(\mu+\cot\phi) x′2(μ+cotϕ).

调整阶数 a a a，使得 μ = − cot ⁡ ϕ \mu=-\cot\phi μ=−cotϕ，消去二次项。

在没有二次项的条件下，再假设 x = f 0 / csc ⁡ ϕ x=f_0/\csc\phi x=f0​/cscϕ，消去一次项，于是积分内部没有 x ′ x' x′，可以去掉积分号，变换结果为 A ϕ exp ⁡ ( i π x 2 cot ⁡ ϕ ) A_\phi\exp(i\pi x^2\cot\phi) Aϕ​exp(iπx2cotϕ)，代入 x = f 0 / csc ⁡ ϕ x=f_0/\csc\phi x=f0​/cscϕ，显然这是一个非零项。

如果 x ≠ f 0 / csc ⁡ ϕ x\neq f_0/\csc\phi x​=f0​/cscϕ，就会出现一个一次项的积分。由于 sin ⁡ t \sin t sint和 cos ⁡ t \cos t cost的周期性，其无穷积分等于零，再根据欧拉公式，有下式成立:

∫ − ∞ ∞ e i t d t = 0 \int^\infty_{-\infty} e^{it}dt=0 ∫−∞∞​eitdt=0

因此一次项的存在会使整个积分为0，换句话说，该积分只有在 x = f 0 / csc ⁡ ϕ x=f_0/\csc\phi x=f0​/cscϕ处有值，而在其它地方为0，类似于冲激函数。

在二次项存在的情况下，该积分就没有上述性质，能量散开，不易检测。

参考文献[1]总结了利用分数阶傅里叶变换检测chirp信号的步骤以及参数估计的公式：

其中 α 0 \alpha_0 α0​在上文中指 ϕ \phi ϕ。整个问题变成了一个二维最优化问题，优化目标是最大化分数阶傅里叶变换函数的值，参数为阶数 a a a和频率 x x x。

注意，在实际操作中发现上述参数估计方程并不正确。进一步查阅文献并参考了他人代码后发现，上述参数估计方程是适用于连续函数的，在实际的计算机中都是离散信号，因此 μ 0 \mu_0 μ0​需要使用归一化调频斜率 k = B f s k=\frac{B}{f_s} k=fs​B​，且 x csc ⁡ ϕ x\csc\phi xcscϕ实际是中心频率而不是起始频率。

仿真

在例如汽车防撞雷达等一些场景中，我们提前知道要检测的chirp信号参数，可以直接根据参数估计分数阶傅里叶分解的阶数，进而直接检测目标信号。

假设我们使用上文中的chirp参数，首先求解归一化调频斜率 k = B f s = 0.0417 k=\frac{B}{f_s}=0.0417 k=fs​B​=0.0417

代入参数估计公式中，计算旋转角度 ϕ = a r c c o t ( − k ) = − 1.5292 \phi=arccot(-k)=-1.5292 ϕ=arccot(−k)=−1.5292，然后计算阶次 a = ϕ / π 2 = − 0.9735 a=\phi/\frac{\pi}{2}=-0.9735 a=ϕ/2π​=−0.9735

负数不影响结果，实际上分数阶傅里叶变换周期为4，而对于chirp信号检测来说， cot ⁡ ( ϕ ) = cot ⁡ ( ϕ + π ) ， \cot(\phi)=\cot(\phi+\pi)， cot(ϕ)=cot(ϕ+π)，周期进一步缩小为2。

对chirp信号做阶数为-0.9735的分数阶傅里叶变换，得到如下时频域图：

由图可以看出，信号在2.8976e10的位置出现了一个非常大的峰值，相比于fft的频域图像，显然这个峰值更容易检测到。因此使用分数阶傅里叶变换比一般的傅里叶变换更抗噪声干扰。

令 x = 2.8976 e 10 x=2.8976e10 x=2.8976e10，估计信号的中心频率为 x csc ⁡ ϕ = 29 G h z x\csc\phi=29Ghz xcscϕ=29Ghz，因此可以推算出起始频率为24Ghz，与预设相符。

代码

close all
fs=24e10;%采样频率
T=1e-7;%时宽
B=10e9;%带宽
mu=B/T;%调频率
n=round(T*fs);%采样点个数
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 on

subplot(212)
plot(t,real(s))%时域图
title("调频信号时域")
xlabel("t/s")
xlim([T-1e-9,T])
grid on

figure
S=fftshift(fft(real(s))./n);
f=linspace(-fs/2,fs/2-1,n);%频域横坐标，注意奈奎斯特采样定理，最大原信号最大频率不超过采样频率的一半
plot(f,abs(S))%频域图
title("发射信号频域")
xlim([-50e9,50e9])
grid on

k=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("发射信号时频域")
% xlim([0,fs/2-1])
title("a="+num2str(a))
grid on
% end

fh_=abs(f(abs(S)==max(abs(S)))*csc(a*pi/2));%估计的中心频率
f0_=fh-B/2;%估计的起始频率
           

参考文献

[1] 渠莹,杨俊. 基于分数阶傅里叶变换的LFM信号参数估计[J]. 物联网技术,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.

注：本人系初学者，以上仅为个人学习笔记，可能存在错误，有任何问题欢迎在评论区讨论。本文将会随着学习理解的加深而继续修改。