Pardon the LaTeX, this is old documentation (and old code). You may be better off just skipping to the references section.
This was a direct port from a c++ implementation (which I've since lost), which is why it's not very 'pythonic'. I'll update it soon, I promise.
import frfft
import numpy as np
t=np.linspace(-4*np.pi,4*np.pi, 1000)
x=np.sin(2*np.pi*40*t) + np.sin(2*np.pi*20*t) + np.sin(2*np.pi*10*t)
X=frfft.FrFFT(x,1./1024)
subplot(211)
plot(np.abs(X.result))
subplot(212)
plot(np.unwrap(arctan2(X.result.imag,X.result.real)))
show()
The FrFFT equation is \f$ G_k(\textbf{\textbf{x}},\alpha) = \displaystyle\sum_{j=0}^{m-1} x_{j}e^{-2 \pi i j k \alpha} \f$.
To calculate it, re-write the equation as
\f$ G_k(\textbf{\textbf{x}}, \alpha) = e^{\pi i k^2 \alpha} \displaystyle\sum_{j=0}^{m-1} y_{j} z_{k-j} \f$
where the m-long sequences y and z are defined by
\f$\y_j = x_j e^{- \pi i j^2 \alpha} \ z_j = e^{\pi i j^2 \alpha} \f$
To ensure a radix-2 FFT, we extend the sequences to a length of \f$ 2p = nextpow2(m+1) \f$ according to the following:
\f$
\begin{tabular}{ l r }
This satisfies the properties for a 2p-point circular convolution and standard FFTs can be used. The formula becomes
\f$ G_k(\textbf{\textit{x}}, \alpha) = e^{- \pi i k^2 \alpha} F_k^{-1} (\textbf{\textit{w}})\f$
Where \f$\textbf{w}\f$ is the 2p-long sequence defined by \f$\textit{w} = F_k(\textbf{y}) F_k(\textbf{z}) \f$
Therefore, the Fractional FFT becomes
\f$ G_k(\textbf{\textit{x}}, \alpha) = e^{- \pi i k^2 \alpha} F_k^{-1} \left(F_k(\textbf{y}) F_k(\textbf{z})\right) \f$
- http://en.wikipedia.org/wiki/Fractional_Fourier_transform
- D. H. Bailey and P. N. Swarztrauber, "The fractional Fourier transform and applications," SIAM Review 33, 389-404 (1991)
- http://www.gnu.org/software/gsl/manual/html_node/Fast-Fourier-Transforms.html