驻定相位原理（POSP）以及线性调频信号的频谱

最近看论文遇到了驻定相位原理，问老师直接给了我一本书让我看，看半天只有一段…不是这个方向的，半路出家做毕业设计需要用到这个定理，有错误的话请不吝赐教。

一、驻定相位原理

在数字信号处理中，经常需要将一个时域信号转换到频域进行处理，比如说需要对其进行傅里叶变换
那么对于这样的一个复数信号：
u(t)=a(t)ejφ(t)u\left(t\right)=a\left(t\right)e^{j\varphi\left(t\right)} u(t)=a(t)ejφ(t)
对其进行傅里叶变换时，就会得到如下式子：
S(ω)=∫−∞∞a(t)ejφ(t)e−jωtdtS(\omega)=\int^\infty_{-\infty}a(t)e^{j\varphi (t)}e^{-j\omega t}dt S(ω)=∫−∞∞a(t)ejφ(t)e−jωtdt
这个积分事实上是很难计算的，并且通常没有解析解。我们现在通常是对时域信号以高于奈奎斯特频率的频率对其进行采样，再进行FFT，利用计算机求得数值解，但在没有计算机的年代，人们只能通过其他方法在特定前提下对其进行估算，比如stationary phase approximation.

This method originates from the 19th century, and is due to George Gabriel Stokes and Lord Kelvin. It is closely related to Laplace’s method and the method of steepest descent, but Laplace’s contribution precedes the others.

对于以下积分式：
P=∫U(t)cosV(t)dt=Re[∫U(t)ejV(t)dt]P=\int U(t)cosV(t)dt=\text{Re} \left[\int U(t)e^{jV(t)}dt\right] P=∫U(t)cosV(t)dt=Re[∫U(t)ejV(t)dt]

给定如下前提，相位V(t)V(t)V(t)对于时间来说是捷变的，幅度U(t)U(t)U(t)是缓变的。对于线性调频信号来说，即
k→∞k\to \infty k→∞

也就是说，当相位随时间变化很快的时候，那么在相位变化一个周期的时间内，幅度可视为一个常数，我们对其积分就近似为零（正负抵消），于是上述积分的值就基本上由V(t)变化缓慢的点决定，即critical points(相位梯度为零的点)。这就是驻定相位原理的基本思想。

二、利用驻定相位原理计算一般调频信号的频谱

设窄带信号的复振幅为
μ(t)=a(t)ejφ(t)\mu(t)=a(t)e^{j\varphi (t)} μ(t)=a(t)ejφ(t)

其傅立叶变换
S(ω)=∫−∞∞a(t)ejφ(t)e−jωtdtS(\omega)=\int^{\infty}_{-\infty}a(t)e^{j\varphi(t)}e^{-j\omega t}dt S(ω)=∫−∞∞a(t)ejφ(t)e−jωtdt

根据驻定相位原理，该积分在
ddt[ωt−φ(t)]=0\frac{d}{dt}[\omega t-\varphi(t)]=0 dtd[ωt−φ(t)]=0

处才显著不为0，用 tkt_ktk 表示驻定相位点（一阶导数为0），则
ω=ϕ′(tk)(∗)\omega=\phi'(t_k) (*) ω=ϕ′(tk)(∗)

对相位项 ωt−ϕ(t)\omega t-\phi(t)ωt−ϕ(t)在驻定相位点附近进行泰勒展开
ωt−φ(t)=ωtk−φ(tk)+[ω−φ′(tk)](t−tk)−φ′′(tk)2(t−tk)2+…\omega t-\varphi (t)=\omega t_k-\varphi(t_k)+\left[\omega-\varphi'(t_k)\right](t-t_k)-\frac{\varphi''(t_k)}{2}(t-t_k)^2+\dots ωt−φ(t)=ωtk−φ(tk)+[ω−φ′(tk)](t−tk)−2φ′′(tk)(t−tk)2+…

忽略二阶以上的高次项（当t远离tkt_ktk时,根据驻定相位原理，积分接近于0，故t必定在tkt_ktk附近取值，即t−tkt-t_kt−tk为小量），并结合（*）式消去一阶项得到：
ωt−φ(t)≈ωtk−φ(tk)−φ′′(tk)2(t−tk)2\omega t-\varphi (t)\approx\omega t_k-\varphi (t_k)-\frac{\varphi''(t_k)}{2}(t-t_k)^2 ωt−φ(t)≈ωtk−φ(tk)−2φ′′(tk)(t−tk)2

代入信号的频域表达式
S(ω)=a(tk)exp{−j[ωtk−φ(tk)]}∫tk−δtk+δexp[jφ′′(tk)2(t−tk)2]dtS(\omega)=a(t_k)\text{exp}\{-j[\omega t_k-\varphi (t_k)]\}\int^{t_k+\delta}_{t_k-\delta}\text{exp}[j\frac{\varphi''(t_k)}{2}(t-t_k)^2]dtS(ω)=a(tk)exp{−j[ωtk−φ(tk)]}∫tk−δtk+δexp[j2φ′′(tk)(t−tk)2]dt

作变量代换
t−tk=u及φ′′(tk)2u2=πy22t-t_k=u及\frac{\varphi''(t_k)}{2}u^2=\frac{\pi y^2}{2} t−tk=u及2φ′′(tk)u2=2πy2

则
du=π[φ′′(tk)]−1/2dydu=\sqrt{\pi}[\varphi''(t_k)]^{-1/2}dy du=π[φ′′(tk)]−1/2dy

代入信号的频域表达式
S(ω)=2πa(tk)φ′′(tk)exp{−j[ωtk−φ(tk)]}∫0φ′′(tk)πδexp(jπy22)dyS(\omega)=2\sqrt{\pi}\frac{a(t_k)}{\sqrt{\varphi''(t_k)}}\text{exp}\{-j[\omega t_k-\varphi (t_k)]\}\int^{\sqrt{\frac{\varphi''(t_k)}{\pi}}\delta}_{0}\text{exp}(j\frac{\pi y^2}{2})dy S(ω)=2πφ′′(tk)a(tk)exp{−j[ωtk−φ(tk)]}∫0πφ′′(tk)δexp(j2πy2)dy

式中积分为菲涅尔积分(Fresnel Integrals可以去查菲涅尔积分表)。若积分上限较大，则菲涅尔积分趋于
12exp(jπ4)\frac{1}{\sqrt{2}}\text{exp}(j\frac{\pi}{4}) 21exp(j4π)

则信号的频域表达式为
S(ω)=2πa(tk)∣φ′′(tk)∣exp[−j(ωtk−φ(tk)−π4)]S(\omega)=\sqrt{2\pi}\frac{a(t_k)}{\sqrt{\left|\varphi''(t_k)\right|}}\text{exp}\left[-j\left(\omega t_k-\varphi(t_k)-\frac{\pi}{4}\right)\right] S(ω)=2π∣φ′′(tk)∣a(tk)exp[−j(ωtk−φ(tk)−4π)]

♣\clubs♣ 驻定相位点不一定只有一个，所以上式是不是应该加个Σ\SigmaΣ?

三、特殊线性调频信号的频谱

所谓线性调频就是

设LFM信号为
u(t)=rect(tTp)exp(jπkrt2)u(t)=rect(\frac{t}{T_p})\text{exp}(j\pi k_r t^2) u(t)=rect(Tpt)exp(jπkrt2)

其傅立叶变换为
U(ω)=∫−∞∞u(t)exp(−jωt)dt=∫−Tp/2Tp/2exp[j(πkrt2−ωt)]dt=exp(−jω24πkr)∫−Tp/2Tp/2exp[j(πkrt−ω2πkr)2]dtU(\omega)=\int^{\infty}_{-\infty}u(t)\text{exp}(-j\omega t)dt =\int^{T_p/2}_{-T_p/2}\text{exp}\left[j\left(\pi k_r t^2-\omega t\right)\right]dt\\ =\text{exp}\left(-j\frac{\omega^2}{4\pi k_r}\right)\int^{T_p/2}_{-T_p/2}\text{exp}\left[j\left(\sqrt{\pi k_r}t-\frac{\omega}{2\sqrt{\pi k_r}}\right)^2\right]dt U(ω)=∫−∞∞u(t)exp(−jωt)dt=∫−Tp/2Tp/2exp[j(πkrt2−ωt)]dt=exp(−j4πkrω2)∫−Tp/2Tp/2exp[j(πkrt−2πkrω)2]dt

作变量代换
(πkrt−ω2πkr)2=π2x2⇒x=2krt−ωπ2krdx=2krdt\left( \sqrt{\pi k_r}t-\frac{\omega}{2\sqrt{\pi k_r}} \right)^2=\frac{\pi}{2}x^2\rArr x=\sqrt{2k_r}t-\frac{\omega}{\pi\sqrt{2k_r}}\\ dx=\sqrt{2k_r}dt (πkrt−2πkrω)2=2πx2⇒x=2krt−π2krωdx=2krdt

则
U(ω)=12krexp(−jω24πkr)∫−X1X2exp(jπx22)dxU(\omega)=\frac{1}{\sqrt{2k_r}}\text{exp}\left( -j\frac{\omega^2}{4\pi k_r} \right) \int^{X_2}_{-X_1}\text{exp} \left( j\frac{\pi x^2}{2} \right)dx U(ω)=2kr1exp(−j4πkrω2)∫−X1X2exp(j2πx2)dx

其中
X1=2krTp2+ωπ2kr=πkrTp+ωπ2krX_1=\sqrt{2k_r}\frac{T_p}{2}+\frac{\omega}{\pi \sqrt{2k_r}}=\frac{\pi k_r T_p+\omega}{\pi \sqrt{2k_r}} X1=2kr2Tp+π2krω=π2krπkrTp+ω
X2=2krTp2−ωπ2kr=πkrTp−ωπ2krX_2=\sqrt{2k_r}\frac{T_p}{2}-\frac{\omega}{\pi \sqrt{2k_r}}=\frac{\pi k_r T_p-\omega}{\pi \sqrt{2k_r}} X2=2kr2Tp−π2krω=π2krπkrTp−ω

菲涅尔积分
C(X)=∫0Xcos⁡(πx22)dxC(X)=\int^{X}_0\cos\left( \frac{\pi x^2}{2}\right)dx C(X)=∫0Xcos(2πx2)dx
S(X)=∫0Xsin⁡(πx22)dxS(X)=\int^{X}_0\sin\left( \frac{\pi x^2}{2}\right)dx S(X)=∫0Xsin(2πx2)dx
C(−X)=−C(X),S(−X)=−S(X)C(-X)=-C(X) , S(-X)= -S(X) C(−X)=−C(X),S(−X)=−S(X)
lim⁡X→∞C(X)=lim⁡X→∞S(X)=0.5\lim_{X\to\infty}C(X)= \lim_{X\to\infty}S(X)=0.5 X→∞limC(X)=X→∞limS(X)=0.5

则矩形包络LFM信号的频谱精确解为
U(ω)=12krexp(−jω24πkr)[C(X1)+jS(X1)+C(X2)+jS(X2)]U(\omega)=\frac{1}{\sqrt{2k_r}}\text{exp}\left( -j\frac{\omega^2}{4\pi k_r}\right) \left[ C(X_1)+jS(X_1)+C(X_2)+jS(X_2) \right] U(ω)=2kr1exp(−j4πkrω2)[C(X1)+jS(X1)+C(X2)+jS(X2)]

其幅度谱和相位谱分别为
∣U(ω)∣=12kr[C(X1)+C(X2)]2+[S(X1)+S(X2)]2\left|U(\omega) \right|=\frac{1}{\sqrt{2k_r}}\sqrt{\left[C(X_1)+C(X_2) \right]^2+\left[S(X_1)+S(X_2) \right]^2} ∣U(ω)∣=2kr1[C(X1)+C(X2)]2+[S(X1)+S(X2)]2
φ(ω)=−ω24πkr+arctan⁡S(X1)+S(X2)C(X1)+C(X2)\varphi(\omega)=-\frac{\omega^2}{4\pi k_r}+\arctan \frac {S(X_1)+S(X_2)} {C(X_1)+C(X_2)} φ(ω)=−4πkrω2+arctanC(X1)+C(X2)S(X1)+S(X2)

在感兴趣的频率范围内，分式
S(X1)+S(X2)C(X1)+C(X2)≈1\frac {S(X_1)+S(X_2)} {C(X_1)+C(X_2)}\approx1 C(X1)+C(X2)S(X1)+S(X2)≈1

当时宽带宽积足够大时，有
∣U(ω)∣=1kr\left|U(\omega) \right|=\frac{1}{\sqrt{k_r}} ∣U(ω)∣=kr1
φ(ω)=−ω24πkr+π4\varphi(\omega)=-\frac{\omega^2}{4\pi k_r}+ \frac{\pi}{4} φ(ω)=−4πkrω2+4π

窄带信号 a(t)e[jϕ(t)]a(t)e^{\left[j\phi(t)\right]}a(t)e[jϕ(t)] 的频谱为
S(ω)=2πa(tk)∣φ′′(tk)∣exp[−j(ωtk−φ(tk)−π4)]S(\omega)=\sqrt{2\pi}\frac{a(t_k)}{\sqrt{\left|\varphi''(t_k)\right|}}\text{exp}\left[-j\left(\omega t_k-\varphi(t_k)-\frac{\pi}{4}\right)\right] S(ω)=2π∣φ′′(tk)∣a(tk)exp[−j(ωtk−φ(tk)−4π)]

对矩形包络的LFM信号（把下面这些特殊信号的取值带到上面的一般结论中）
a(t)=rect(tTp)φ(t)=πkrt2⇒φ′(t)=2πkrt⇒φ′′(t)=2πkrddt[φ(t)−ωt]=0⇒2πkrt−ω=0⇒tk=ω2πkra(t)=rect(\frac{t}{T_p})\\ \varphi(t)=\pi k_r t^2 \rArr \varphi'(t)=2\pi k_r t\rArr \varphi''(t)=2\pi k_r\\ \frac{d}{dt}\left[ \varphi(t)-\omega t\right]=0\rArr 2\pi k_r t-\omega=0\rArr t_k=\frac{\omega}{2\pi k_r} a(t)=rect(Tpt)φ(t)=πkrt2⇒φ′(t)=2πkrt⇒φ′′(t)=2πkrdtd[φ(t)−ωt]=0⇒2πkrt−ω=0⇒tk=2πkrω

利用驻定相位原理得到LFM信号的幅度和相位谱分别为
∣S(ω)∣=2πrect(ω2πkrTp)2πkr=1krrect(ωΔω)\left| S(\omega) \right|=\sqrt{2\pi}\frac{rect\left(\frac{\omega}{2\pi k_r T_p}\right)}{\sqrt{2\pi k_r}}=\frac{1}{\sqrt{k_r}}rect(\frac{\omega}{\Delta \omega})\\ ∣S(ω)∣=2π2πkrrect(2πkrTpω)=kr1rect(Δωω)
Φ(ω)=−ωtk+φ(tk)+π4=−ωω2πkr+πkr(ω2πkr)2+π4=−ω24πkr+π4\Phi(\omega)=-\omega t_k+\varphi(t_k)+\frac{\pi}{4}=-\omega\frac{\omega}{2\pi k_r}+\pi k_r \left( \frac{\omega}{2\pi k_r} \right)^2+\frac{\pi}{4}=-\frac{\omega^2}{4\pi k_r}+\frac{\pi}{4} Φ(ω)=−ωtk+φ(tk)+4π=−ω2πkrω+πkr(2πkrω)2+4π=−4πkrω2+4π

四、快速POSP

设g(t)g(t)g(t)是一个调频信号
g(t)=w(t)exp[jϕ(t)]g(t)=w(t)\mathrm{exp}[j\phi(t)] g(t)=w(t)exp[jϕ(t)]
其中w(t)w(t)w(t)是实包络，ϕ(t)\phi(t)ϕ(t)为信号调制相位。假设与相位相比，包络为时间缓变函数。
该信号的频谱为
G(f)=∫−∞+∞g(t)exp[−j2πft]dt=∫−∞+∞w(t)exp[jϕ(t)−2πft]dt=∫−∞+∞w(t)exp[jθ(t)]dt\begin{aligned} G(f)&=\int_{-\infty}^{+\infty}g(t)\mathrm{exp}[-j2\pi ft]dt\\&=\int_{-\infty}^{+\infty}w(t)\mathrm{exp}[j\phi(t)-2\pi ft]dt\\&=\int_{-\infty}^{+\infty}w(t)\mathrm{exp}[j\theta(t)]dt \end{aligned} G(f)=∫−∞+∞g(t)exp[−j2πft]dt=∫−∞+∞w(t)exp[jϕ(t)−2πft]dt=∫−∞+∞w(t)exp[jθ(t)]dt
最终，G(f)G(f)G(f)的频谱形式如下：
G(f)≈C1W(f)exp[j(Θ(f)±π4)]G(f)\approx C_1W(f)\mathrm{exp}[j(\Theta (f)\pm\frac{\pi}{4})] G(f)≈C1W(f)exp[j(Θ(f)±4π)]
其中变量定义如下：
♣\color{salmon}\clubsuit♣ C1C_1C1为一个通常可忽略的常数
C1=2π∣ϕ′′(t)∣C_1=\sqrt{\frac{2\pi}{\lvert\phi^{''}(t)\rvert}} C1=∣ϕ′′(t)∣2π
♣\color{salmon}\clubsuit♣ W(f)W(f)W(f)为频域包络
W(f)=w[t(f)]W(f)=w[t(f)] W(f)=w[t(f)]
♣\color{salmon}\clubsuit♣ Θ(f)\Theta(f)Θ(f)为频域相位
Θ(f)=θ[t(f)]\Theta(f)=\theta[t(f)] Θ(f)=θ[t(f)]
♣\color{salmon}\clubsuit♣ t(f)t(f)t(f)由信号时频关系给出
dθ(t)dt=0\frac{d\theta(t)}{dt}=0 dtdθ(t)=0
♣\color{salmon}\clubsuit♣ π4\frac{\pi}{4}4π的符号由ϕ′′(t)\phi^{''}(t)ϕ′′(t)的符号给出。与C1C_1C1类似，在绝大多数分析中都可忽略该常数相位的影响。

五、利用快速POSP求解LFM的频谱

LFM信号频谱为
G(f)=∫−∞∞rect(tT)exp(jπKt2)exp(−j2πft)dtG(f)=\int^{\infty}_{-\infty}rect(\frac{t}{T})\text{exp}(j\pi K t^2)\text{exp}(-j2\pi f t)dt G(f)=∫−∞∞rect(Tt)exp(jπKt2)exp(−j2πft)dt
其中包络w(t)w(t)w(t)为矩形，被积相位为
θ(t)=πKt2−2πft\theta(t)=\pi Kt^2-2\pi f t θ(t)=πKt2−2πft
对其求导，并令为0，可得时频关系：
t=fKt=\frac{f}{K} t=Kf
频域相位为
Θ(f)=θ(t=f/K)=πK(fK)2−2πf(fK)=−πf2K\Theta(f)=\theta(t=f/K)=\pi K(\frac{f}{K})^2-2\pi f(\frac{f}{K})=-\pi\frac{f^2}{K} Θ(f)=θ(t=f/K)=πK(Kf)2−2πf(Kf)=−πKf2
频域包络为
W(f)=w(t=f/K)=rect(f∣K∣T)W(f)=w(t=f/K)=\mathrm{rect}(\frac{f}{\lvert K\rvert T}) W(f)=w(t=f/K)=rect(∣K∣Tf)
忽略前面的系数和线性相位，可得频谱为：
G(f)=rect(fKT)exp(−jπf2K)G(f)=\mathrm{rect}(\frac{f}{KT})\mathrm{exp}(-j\pi\frac{f^2}{K}) G(f)=rect(KTf)exp(−jπKf2)
可以看到与前面的结论是一致的。

附录（我的草稿）

参考文献

（1）https://wenku.baidu.com/view/0f5ef8f74693daef5ef73dfd.html
（2）https://en.wikipedia.org/wiki/Chirp_spectrum
（3）https://en.wikipedia.org/wiki/Stationary_phase_approximation
（4）合成孔径雷达成像算法与实现/IAN G.CUMMING