The Fourier transform of sinusoids with cubic phase evolution


If you Fourier transform a quasi-sinusoidal signal with a quadratic phase evolution then the Fourier transform can be approximated using the Fresnel integrals. However, what if your signal has a cubic phase evolution? Can you go beyond the Fresnel integrals? This is to expand on the work of Davies (2015).

In this particular case we are interested in signal that may have intrinsic cubic phase evolution, but the majority of any cubic component is from Doppler shifting of the signal e.g. in a binary system orbit.

Firstly, we will define the phase evolution of our signal as

\begin{equation} \phi(t)\approx 2\pi\sum_{l=0}\frac{f^{(l)}}{(l+1)!}\left(t+\delta(t)-t_{0}\right)^{(l+1)},\\ \end{equation}

where \(f^{(l)}\) are the \(l^{\rm th}\) intrinsic frequency derivatives of the signal defined at the epoch \(t_{0}\), and \(\delta(t)\) is the time delay between the source inertial frame and the solar system barycentre (SSB).

If we Taylor expand \(\delta(t)\) to second order we get

\begin{equation} \delta(t)\approx\delta_{k}+\dot{\delta_{k}}t+\frac{1}{2}\ddot{\delta_{k}}t^{2},\\ \end{equation}

where \(\delta_{k}=\delta(t_{k})\) is the time delay at time \(t_{k}\), \(\dot{\delta_{k}}={\rm d}\delta(t_{k})/{\rm d}t\) and \(\ddot{\delta_{k}}={\rm d}^{2}\delta(t_{k})/{\rm d}t^{2}\).

Therefore we can expand out the phase evolution to be

\begin{equation} \label{eq:phaseevo}\phi(t)\approx 2\pi\sum_{l=0}\frac{f^{(l)}}{(l+1)!}\left(t\left(1+\dot{\delta_{k}}+\frac{\ddot{\delta_{k}}}{2}t\right)-t_{0}+\delta_{k}+t_{k}\right)^{l+1}.\\ \end{equation}

We want to have equation \ref{eq:phaseevo} in the following form

\begin{equation} \phi(t)=\phi_{k}+2\pi\left(f_{k}t+\frac{\dot{f_{k}}}{2}t^{2}+\frac{\ddot{f_{k}}}{6}t^{3}\right),\\ \end{equation}

so to find \(\phi_{k}\), \(f_{k}\), \(\dot{f_{k}}\) and \(\ddot{f_{k}}\) we need to Taylor expand equation \ref{eq:phaseevo} about \(t=0\). This gives the following values

\begin{align} \phi_{k} & =2\pi\sum_{l=0}\frac{f^{(l)}}{(l+1)!}\left(t_{k}-t_{0}+\delta_{k}\right)^{l+1} \\ f_{k} & =(1+\dot{\delta_{k}})\sum_{l=0}\frac{f^{(l)}}{l!}\left(t_{k}-t_{0}+\delta_{k}\right)^{l} \\ \dot{f_{k}} & =\sum_{l=0}\frac{f^{(l)}}{l!}\left(\ddot{\delta_{k}}+\frac{l(1+\dot{\delta_{k}})^{2}}{t_{k}-t_{0}+\delta_{k}}\right)\left(t_{k}-t_{0}+\delta_{k}\right)^{l} \\ & =\ddot{\delta_{k}}\left(\sum_{l=0}\frac{f^{(l)}}{l!}\left(t_{k}-t_{0}+\delta_{k}\right)^{l}\right)+(1+\dot{\delta_{k}})^{2}\left(\sum_{l=0}\frac{f^{(l)}}{(l-1)!}\left(t_{k}-t_{0}+\delta_{k}\right)^{l-1}\right)\\ \end{align}

whilst \(\ddot{f_{k}}\) can be approximated as

\begin{equation} \ddot{f_{k}}=\frac{\dot{f}_{\rm end}-\dot{f}_{\rm start}}{\Delta t},\\ \end{equation}

where \(t_{k}\) is the central time of the Fourier transformed data, \(\Delta t\) is the length of the data, and

\begin{align} \dot{f}_{\rm start/end}= & \ddot{\delta}(t_{\rm s/e})\left(\sum_{l=0}\frac{f^{(l)}}{l!}\left(t_{\rm s/e}-t_{0}+\delta(t_{\rm s/e})\right)^{l}\right)\notag \\ & +(1+\dot{\delta}(t_{\rm s/e}))^{2}\left(\sum_{l=0}\frac{f^{(l)}}{(l-1)!}\left(t_{\rm s/e}-t_{0}+\delta(t_{\rm s/e})\right)^{l-1}\right),\\ \end{align}

where \(t_{\rm s}=t_{k}-\Delta t/2\) and \(t_{\rm e}=t_{k}+\Delta t/2\).

We can define the Fourier transform of our signal, which has complex amplitude at time \(t_{k}\) of \(y_{k}\) as

\begin{align} H_{k}(f)\approx & y_{k}e^{i\phi_{k}-i\pi f\Delta t}\int_{-\Delta t/2}^{\Delta t/2}\exp{\left[2\pi i\left((f_{k}-f)t+\frac{t^{2}}{2}\dot{f_{k}}+\frac{t^{3}}{6}\ddot{f_{k}}\right)\right]}{\rm d}t\notag \\ & \label{eq:fourier}+y_{k}^{*}e^{-i\phi_{k}-i\pi f\Delta t}\int_{-\Delta t/2}^{\Delta t/2}\exp{\left[-2\pi i\left((f_{k}+f)t-\frac{t^{2}}{2}\dot{f_{k}}-\frac{t^{3}}{6}\ddot{f_{k}}\right)\right]}{\rm d}t.\\ \end{align}

The second term in this equation involving \((f+f_{k})\) is negligible for our purposes, so only the first term is required.

We want to convert this into a more easily solvable form e.g.

\begin{equation} \label{eq:fourier2}H_{k}(f)=y_{k}^{\prime}e^{i\phi_{k}-i\pi f\Delta t}\int_{z_{1}}^{z_{2}}\exp{\left[i\left(z^{3}+\alpha z+\beta\right)\right]}{\rm d}z.\\ \end{equation}

Firstly, lets put the cubic in equation \ref{eq:fourier} in the form

\begin{equation} \label{eq:subst}M=at^{3}+bt^{2}+ct+d,\\ \end{equation}

where \(a=(\pi/3)\ddot{f_{k}}\), \(b=\pi\dot{f_{k}}\), \(c=2\pi(f_{k}-f)\) and \(d=0\) and make the substitution \(t=Az+B\). To get the derivative \({\rm d}z\) we have

\begin{equation} \frac{{\rm d}t}{{\rm d}z}=A,\\ \end{equation}

so, \({\rm d}t=A{\rm d}z\) and therefore in equation \ref{eq:fourier2} we have \(y_{k}^{\prime}=Ay_{k}\). If we make the substitution into equation \ref{eq:subst} then we get

\begin{equation} \label{eq:Mexpanded}M=A^{3}az^{3}+A^{2}(3Ba+b)z^{2}+A(3B^{2}a+2Bb+c)z+B(B^{2}a+Bb+c)+d.\\ \end{equation}

If we equation the terms in equation \ref{eq:Mexpanded} with those in exponential term within the integral in equation \ref{eq:fourier2} we get

\begin{align} A^{3}a= & 1,\notag \\ A^{2}(3Ba+b)= & \label{eq:firstterms} 0,\\ \end{align}


\begin{align} \alpha= & A(3B^{2}a+2Bb+c),\notag \\ \beta= & \label{eq:secondterms} B(B^{2}a+Bb+c)+d.\\ \end{align}

We can use equation \ref{eq:firstterms} to solve for \(A\) and \(B\), with the real roots giving

\begin{align} A= & a^{-1/3},\notag \\ B= & -\frac{b}{3a}.\\ \end{align}

Substituting for the values of \(a\) and \(b\) we have

\begin{equation} \label{eq:A}A=\left(\frac{\pi}{3}\ddot{f_{k}}\right)^{-1/3}\\ \end{equation}


\begin{equation} \label{eq:B}B=-\left(\frac{\pi}{3}\right)^{4/3}\dot{f_{k}}\ddot{f_{k}}^{1/3}.\\ \end{equation}