基于分数傅里叶变换的chirp信号参数估计-易微帮

系列文章目录

《SAR信号处理重要工具-傅里叶变换》

《分数傅里叶变换》

文章目录

前言

《SAR信号处理重要工具-傅里叶变换》较为详细的介绍了傅里叶变换的来源以及一些使用技巧。《分数傅里叶变换》重点介绍了信号领域比较常见的一种信号类型——chirp信号，并利用FFT简单分析了chirp信号的频域特征，在时宽带宽积较大时，其频谱近似矩形。通过对chirp信号进行短时傅里叶变换，不难发现chirp信号频率时变特性。为了更好了解chirp信号的特点，我们利用分数阶傅里叶变换分析了chirp信号，发现其在特定阶次上表现出很好的聚焦特性。基于这一特性，本文简单介绍基于分数傅里叶变换的chirp信号参数估计。

一、参数估计原理

分数傅里叶变换的一种表达式为：

$X_{\alpha }\left ( u \right )=A_{\alpha }e^{j\pi \cot\alpha u^{2}}\int_{-\infty }^{\infty }x\left ( t \right )e^{j\pi\cot\alpha t^{2}}e^{-j2\pi\csc\alpha tu}dt$

其中 $A_{\alpha }=\sqrt{\frac{1-j\cot\alpha }{2\pi}}$ ， $\alpha=p\frac{\pi}{2}$ ， $p$ 表示阶次， $u$ 表示分数阶傅里叶域，当阶次为1时，上式退化为傅里叶变换。

当信号 $x\left ( t \right )$ 为：

$x\left ( t \right )=Ae^{j\pi Kt^{2}+j2\pi ft+\phi },-T/2\leq t\leq T/2$

则：

$X_{\alpha }\left ( u \right )=A_{\alpha }e^{j\pi \cot\alpha u^{2}}Ae^{j\phi }\int_{-T/2 }^{T/2 }e^{j\pi\left ( \cot\alpha+K \right ) t^{2}}e^{j2\pi t\left ( f-u\csc\alpha \right )}dt$

不难发现，当 $K+\cot\alpha =0,f-u\csc\alpha =0$ 时

$\max\left \{ \left |X_{\alpha }\left ( u \right ) \right | \right \}=\left |A_{\alpha } \right |AT$

因此可以根据搜索变量 $\left (\alpha, u \right )$ 下搜索函数 $\left |X_{\alpha }\left ( u \right ) \right |$ 的峰值位置 $\left (\hat{\alpha}, \hat{u} \right )$ 以及该位置处的变换结果 $X_{\hat{\alpha} }\left ( \hat{u} \right )$ 得到chirp信号幅度，初相，频率以及调频类的估计结果：

$\hat{A}=\frac{\left | X_{\hat{\alpha} }\left ( \hat{u} \right )\right |}{T\left | A_{\hat{\alpha }} \right |}$ ， $\hat{\phi }=\arg\left [\frac{X_{\hat{\alpha} }\left ( \hat{u} \right )}{A_{\hat{\alpha }} e^{j\pi \cot\hat{\alpha} { \hat{u}}^{2}}} \right ]$ ， $\hat{f}=\hat{u} \csc\hat{\alpha }$ ， $\hat{K}=- \cot\hat{\alpha }$

二、离散分数阶傅里叶变换

考虑到实际处理的信号一般为数字信号，为了将分数阶傅里叶变换更好的应用于数字信号中，这一节内容将详细介绍离散分数阶傅里叶变化。在介绍分数傅里叶变换的离散算法之前，先来介绍量纲归一化。

2.1、量纲归一化

这里对信号进行量纲归一化的目的主要有两个：首先通过量纲归一化，使得各分数域支撑域一致，便于二维搜索结果的展示；此外各分数域统一，便于之后构造离散算法的卷积形式，从而借助傅里叶变换的快速算法实现对分数傅里叶的快速计算。上图显示了量纲归一化的原理，从右图可以看出量纲归一化后，信号的时宽、带宽一致，各分数域的支撑域一致。归一化之前离散信号时间宽度 $T$ s，频域宽度为采样率 $f_{s}$ Hz，量纲归一化（尺度和单位同时归一化）将信号时长伸缩到 $\Delta x$ ，采样率变为 $\Delta x$ 。

这种归一化的操作虽然使得各分数域的支撑域一致，但是，也破坏了原信号的结构，以chirp信号为例，由于时频轴放缩的原因，归一化后的调频率变成了原信号调频率的S平方倍，频率变为原频率的S倍。因此量纲归一化后的信号进行离散分数阶傅里叶变换（基于FFT的快速算法）估计参数时，需要考虑这方面的影响，即最终估计参数表达式为：

$\hat{A}=\frac{\left | X_{\hat{\alpha} }\left ( \hat{u} \right )\right |}{T\left | A_{\hat{\alpha }} \right |}$ ， $\hat{\phi }=\arg\left [\frac{X_{\hat{\alpha} }\left ( \hat{u} \right )}{A_{\hat{\alpha }} e^{j\pi \cot\hat{\alpha} { \hat{u}}^{2}}} \right ]$ ， $\hat{f}=\frac{\hat{u} \csc\hat{\alpha }}{S}$ ， $\hat{K}=- \frac{\cot\hat{\alpha }}{S^{2}}$

仔细观察不难发现，量纲归一化是对模拟信号的影响，量纲归一化前后的离散信号并没有区别，并且这一结论可以由下面的仿真进一步验证。

2.2、离散分数阶傅里叶变换的快速实现

量纲归一化的一个目的是构造离散算法的卷积形式，从而借助傅里叶变换的快速算法实现对分数傅里叶的快速计算。对量纲归一化的信号 $x\left ( u_{0} \right )$ 进行分数阶傅里叶变换：

$x_{\alpha }\left ( u_{\alpha } \right )=A_{\alpha }e^{j\pi \cot\alpha u_{\alpha }^{2}}\int_{-\infty }^{\infty }x\left ( u_{0} \right )e^{j\pi\cot\alpha u_{0}^{2}}e^{-j2\pi\csc\alpha u_{0}u_{\alpha }}du_{0}$

考虑到实际的信号是由模拟信号采样得到的数字信号，因此可以由香农插值定理，将 $x\left ( u_{0} \right )e^{j\pi\cot\alpha u_{0}^{2}}$ 用对应的采样值表示，考虑到 $x\left ( u_{0} \right )$ 和 $e^{j\pi\cot\alpha u_{0}^{2}}$ 在 $\left |\cot\alpha \right |\leq 1$ 情况下带宽都小于或等于 $\Delta x$ ，根据信号在时域上相乘等效为信号在频域上卷积后IFFT的结果，可以判断 $x\left ( u_{0} \right )e^{j\pi\cot\alpha u_{0}^{2}}$ 的带宽不大于 $2\Delta x$ ，因此 $x\left ( u_{0} \right )e^{j\pi\cot\alpha u_{0}^{2}}$ 的采样率需要为 $2\Delta x$ 才能由离散信号重构出离散信号。而 $x\left ( u_{0} \right )$ 只有 $\Delta x$ 采样率的离散信号 $x\left ( \frac{n}{\Delta x} \right ),n=0,1,\cdots ,N-1$ ，可以对采样率为 $\Delta x$ 的采样信号以采样率 $2\Delta x$ 进行间隔补零，然后在频域上进行 $\Delta x$ 带宽的低通滤波处理得到采样率为 $2\Delta x$ 的采样信号 $x\left ( \frac{n}{2\Delta x} \right ),n=0,1,\cdots ,2N-1$ 。

$x\left ( u_{0} \right )e^{j\pi\cot\alpha u_{0}^{2}}=\sum_{n=0}^{2N-1}x\left ( \frac{n}{2\Delta x} \right )e^{j\pi \cot\alpha \left ( \frac{n}{2\Delta x} \right )^{2}}sinc\left [ 2\Delta x\left ( u_{0}- \frac{n}{2\Delta x}\right ) \right ]$

将重构的信号代入分数傅里叶变换公式得：

$x_{\alpha }\left ( u_{\alpha } \right )=A_{\alpha }e^{j\pi \cot\alpha u_{\alpha }^{2}}\int_{-\infty }^{\infty }\sum_{n=0}^{2N-1}x\left ( \frac{n}{2\Delta x} \right )e^{j\pi \cot\alpha \left ( \frac{n}{2\Delta x} \right )^{2}}sinc\left [ 2\Delta x\left ( u_{0}- \frac{n}{2\Delta x}\right ) \right ]e^{-j2\pi\csc\alpha u_{0}u_{\alpha }}du_{0}$

可以看出变换公式中同时出现积分、求和运算，这里我们交换积分、求和的运算顺序，并且这种交换不会影响最后结果的计算：

$x_{\alpha }\left ( u_{\alpha } \right )=A_{\alpha }e^{j\pi \cot\alpha u_{\alpha }^{2}}\sum_{n=0}^{2N-1} \int_{-\infty }^{\infty }sinc\left [ 2\Delta x\left ( u_{0}- \frac{n}{2\Delta x}\right ) \right ]e^{-j2\pi\csc\alpha u_{0}u_{\alpha }}du_{0}x\left ( \frac{n}{2\Delta x} \right )e^{j\pi \cot\alpha \left ( \frac{n}{2\Delta x} \right )^{2}}$

积分部分是时移的sinc函数的伸缩傅里叶变换，通过变量替换以及傅里叶变换规律可以很快得到积分结果：

$\int_{-\infty }^{\infty }sinc\left [ 2\Delta x\left ( u_{0}- \frac{n}{2\Delta x}\right ) \right ]e^{-j2\pi\csc\alpha u_{0}u_{\alpha }}du_{0}=\frac{1}{2\Delta x} rect\left [ \frac{u_{\alpha }\csc\alpha }{2\Delta x} \right ] e^{-j2\pi\csc\alpha u_{\alpha }\frac{n}{2\Delta x}}$

由此

$x_{\alpha }\left ( u_{\alpha } \right )=\frac{A_{\alpha }}{2\Delta x} e^{j\pi \cot\alpha u_{\alpha }^{2}}\sum_{n=0}^{2N-1}x\left ( \frac{n}{2\Delta x} \right )e^{j\pi \cot\alpha \left ( \frac{n}{2\Delta x} \right )^{2}}e^{-j2\pi\csc\alpha u_{\alpha }\frac{n}{2\Delta x}}$

对变换结果在该分数域上进行离散取值，离散的间隔为2*delta分之一:

$x_{\alpha }\left ( \frac{m}{2\Delta x} \right )=\frac{A_{\alpha }}{2\Delta x} \sum_{n=0}^{2N-1}x\left ( \frac{n}{2\Delta x} \right )e^{j\pi \cot\alpha \left ( \frac{n}{2\Delta x} \right )^{2}+j\pi \cot\alpha \left ( \frac{m}{2\Delta x} \right )^{2}-j2\pi\csc\alpha \frac{mn}{\left (2\Delta x \right )^{2}}}$

构造离散信号的卷积形式：

$x_{\alpha }\left ( \frac{m}{2\Delta x} \right )=\frac{A_{\alpha }}{2\Delta x} e^{j\pi\left ( \cot\alpha -\csc\alpha \right ) \left ( \frac{m}{2\Delta x} \right )^{2}}\sum_{n=0}^{2N-1}x\left ( \frac{n}{2\Delta x} \right )e^{j\pi\left ( \cot\alpha -\csc\alpha \right ) \left ( \frac{n}{2\Delta x} \right )^{2}}e^{j\pi\csc\alpha \frac{\left (m-n \right )^{2}}{\left (2\Delta x \right )^{2}}}$

由此可借助傅里叶变换的快速算法实现对分数傅里叶变换的快速计算。

三、仿真分析

3.1、单分量

单个chirp信号，调频类200 Hz/s，频率100 Hz，初相相位1 rad，幅度10，信号时长1 s，采样率1000 Hz。

为了得到调频率的高精度估计结果，同时考虑计算效率的问题，仿真对分数傅里叶变换的阶次采取先粗搜后精搜的策略，从而实现对阶次的高精度估计，之后根据调频率、频率与搜索峰值位置之间的关系得到调频率与频率的估值。最后调频率、频率的估值为：

$\hat{K}=\frac{-\cot \hat{\alpha }}{S^{2}}\approx 200.3824 \text{Hz/s}$

$\hat{f}=\frac{\hat{u}\csc\hat{\alpha }}{S}\approx 99.9481 \text{Hz}$

3.2、多分量

两个chirp信号叠加：调频类200 Hz/s，频率100 Hz，初相相位1 rad，幅度10；调频类400 Hz/s，频率200 Hz，初相相位1 rad，幅度10。信号时长1 s，采样率1000 Hz。

通过粗搜确定原始信号含有两个chirp信号分量，并得到两个chirp信号的大致位置。分别在这两个chirp信号分量对应的阶次附近进行精搜。最后根据调频率、频率与搜索峰值位置之间的关系得到这两个chirp信号分量的调频率与频率的估值:

$\hat{K_{1}}=\frac{-\cot \hat{\alpha_{1} }}{S^{2}}\approx 200.3824 \text{Hz/s}\\ \hat{K_{2}}=\frac{-\cot \hat{\alpha_{2} }}{S^{2}}\approx 399.9310 \text{Hz/s}$

$\hat{f_{1}}=\frac{\hat{u_{1}}\csc\hat{\alpha_{1} }}{S}\approx 99.9481 \text{Hz}\\ \hat{f_{2}}=\frac{\hat{u_{2}}\csc\hat{\alpha_{2} }}{S}\approx 200.3234\text{Hz}$