# Debug the difference between two Frft implementations

I need help to understand the difference between two fractional Fourier transform implementations.

On this website two different implementations for the fractional Fourier transform are presented (frft, and frft2). The two implementations are provided as Matlab code.

The accompanying paper on the linked website describes the second implementation as using less resources but delivering an identical result.

I want to visualize the Fractional Fourier transform inside a web page. Of course I would prefer to use the more performant implementation. This is why I tried to translate the second implementation (frft2) into JavaScript.

While doing so I noticed that the resulting signal is almost but not exactly what I expected. For example I expected an even rect signal to stay approximately even when being frft-transformed (as described on wikipedia).

To double check I also translated the first (frft) implementation and noticed that now the results are as expected. So either the two implementations differ in their output, or I made a mistake translating the one implementation but not the other.

To double check again and also to stay closer to the syntax and semantics of the Matlab code, I also translated both implementations into Python using NumPy. Unfortunately this resulted in the same difference: frft working as expected and frft2 being slightly off (even more off for time shifted signals).

Below you can see a plot comparing real part of frft(rect(t), a=0.49) (green) and frft2(rect(t), a=0.49) (red):

In this plot the difference is not that big, but depending on the time shift of the rect and on the fractional parameter a, both implementations differ quite a lot.

This is a Jupyter notebook containing both my Python implementations.

This is my Python implementation of frft working correctly (I think):

# Matlab to Python translation from
# https://nalag.cs.kuleuven.be/research/software/FRFT/frft.m

import numpy as np
import scipy
import scipy.signal

def frft(f, a):
ret = np.zeros_like(f, dtype=np.complex128)
f = f.copy().astype(np.complex128)
N = len(f)
shft = np.fmod(np.arange(N) + np.fix(N / 2), N).astype(int)
sN = np.sqrt(N)
a = np.remainder(a, 4.0)

if a == 0.0:
return f
if a == 2.0:
return np.flipud(f)
if a == 1.0:
ret[shft] = np.fft.fft(f[shft]) / sN
return ret
if a == 3.0:
ret[shft] = np.fft.ifft(f[shft]) * sN
return ret

# reduce to interval 0.5 < a < 1.5
if a > 2.0:
a = a - 2.0
f = np.flipud(f)
if a > 1.5:
a = a - 1
f[shft] = np.fft.fft(f[shft]) / sN
if a < 0.5:
a = a + 1
f[shft] = np.fft.ifft(f[shft]) * sN

# the general case for 0.5 < a < 1.5
alpha = a * np.pi / 2
tana2 = np.tan(alpha / 2)
sina = np.sin(alpha)
f = np.hstack((np.zeros(N - 1), sincinterp(f), np.zeros(N - 1))).T

# chirp premultiplication
chrp = np.exp(-1j * np.pi / N * tana2 / 4 * np.arange(-2 * N + 2, 2 * N - 1).T ** 2)
f = chrp * f

# chirp convolution
c = np.pi / N / sina / 4
ret = scipy.signal.fftconvolve(
np.exp(1j * c * np.arange(-(4 * N - 4), 4 * N - 3).T ** 2), f
)
ret = ret[4 * N - 4 : 8 * N - 7] * np.sqrt(c / np.pi)

# chirp post multiplication
ret = chrp * ret

# normalizing constant
ret = np.exp(-1j * (1 - a) * np.pi / 4) * ret[N - 1 : -N + 1 : 2]

return ret

def ifrft(f, a):
return frft(f, -a)

def sincinterp(x):
N = len(x)
y = np.zeros(2 * N - 1, dtype=x.dtype)
y[: 2 * N : 2] = x
xint = scipy.signal.fftconvolve(
y[: 2 * N],
np.sinc(np.arange(-(2 * N - 3), (2 * N - 2)).T / 2),
)

return xint[2 * N - 3 : -2 * N + 3]


This is my Python implementation of frft2 being slightly off:

# Matlab to Python translation from
# https://nalag.cs.kuleuven.be/research/software/FRFT/frft2.m

import numpy as np
import scipy
import scipy.signal

def frft2(f, a):
f0 = f.flatten()
N = len(f)
sN = np.sqrt(N)
a = np.mod(a, 4)

if a == 0:
return f0
elif a == 2:
return np.flipud(f0)
elif a == 1:
return np.fft.ifftshift(np.fft.fft(np.fft.fftshift(f0))) / sN
elif a == 3:
return np.fft.ifftshift(np.fft.ifft(np.fft.fftshift(f0))) * sN

if a > 2.0:
a = a - 2
f0 = np.flipud(f0)
if a > 1.5:
a = a - 1
f0 = np.fft.ifftshift(np.fft.fft(np.fft.fftshift(f0))) / sN
if a < 0.5:
a = a + 1
f0 = np.fft.ifftshift(np.fft.ifft(np.fft.fftshift(f0))) * sN

alpha = a * np.pi / 2
s = np.pi / (N + 1) / np.sin(alpha) / 4
t = np.pi / (N + 1) * np.tan(alpha / 2) / 4
Cs = np.sqrt(s / np.pi) * np.exp(-1j * (1 - a) * np.pi / 4)

f1 = fconv(f0, np.sinc(np.arange(-(2 * N - 3), 2 * N - 1, 2) / 2), 1)
f1 = f1[N : 2 * N]
chrp = np.exp(-1j * t * np.arange(-N, N) ** 2)
l0 = chrp[::2]
l1 = chrp[1::2]
f0 = f0 * l0
f1 = f1 * l1
chrp = np.exp(1j * s * np.arange(-(2 * N), 2 * N - 1) ** 2)
e1 = chrp[::2]
e0 = chrp[1::2]
f0 = fconv(f0, e0, 0)
f1 = fconv(f1, e1, 0)
h0 = np.fft.ifft(f0 + f1)

Faf = Cs * l0 * h0[N : 2 * N]

return Faf

def fconv(x, y, c):
N = len(np.concatenate((x.flatten(), y.flatten()))) - 1
P = 2 ** np.ceil(np.log2(N)).astype(int)
z = np.fft.fft(x, P) * np.fft.fft(y, P)

if c != 0:
z = np.fft.ifft(z)
z = z[N::-1]

return z


Since I have not much experience in Matlab I suspect that I mistranslated some of the range syntax and included some off-by-1 errors. I already spent some hours debugging.

Both implementations map the parameter a into the range (0.5, 1.5) to make use of symmetries but for the frft2 implementation I noticed a discontinuity at a=0.5 and at a=2.5. To me it looks like this case might be handle incorrectly:

if a < 0.5:
a = a + 1
f0 = np.fft.ifftshift(np.fft.ifft(np.fft.fftshift(f0))) * sN


For example with frft an even real signal is even for both a=0.45 and a=0.55. But with frft2 an even real signal is complex and odd for 0<a<=0.5 but suddenly even for a>0.5.

But I could not find the true reason for this and since both implementations map a into this range the actual difference might be somethere else.

I would be glad if someone could either point out my mistake or confirm that the two algorithms actually differ in their results.

Update: Now I also tried to run the Matlab code (in Octave). It looks like even the original Matlab code leads to different results. But I still do not understand why. Numerical instability? Is one of the algorithms simply not correct? I am mostly confused by the difference in aproximate odd/eveness of the results and the incontinuities regarding a.

octave> frft([1,0,0,0], 0.50)
0.003910 - 0.548475i
0.393301 - 0.072756i
-0.349355 + 0.152855i
0.051604 + 0.053868i

octave> frft2([1,0,0,0], 0.50)
-0.0753 - 0.4989i
0.2860 + 0.1183i
-0.1111 + 0.1219i
0.3095 - 0.1835i

octave> frft([1,0,0,0], 0.51)
-0.012891 - 0.545906i
0.390336 - 0.081043i
-0.340700 + 0.175301i
0.081861 + 0.048968i

octave> frft2([1,0,0,0], 0.51)
-0.0962 - 0.4873i
0.2876 + 0.1104i
-0.0996 + 0.1328i
0.3034 - 0.2179i

octave> frft([1,0,0,0], 0.49)
0.3110 - 0.2646i
-0.2468 + 0.2129i
0.2468 - 0.2129i
-0.3110 + 0.2646i

octave> frft2([1,0,0,0], 0.49)
0.3207 - 0.0250i
-0.2658 - 0.0597i
0.2658 + 0.0597i
-0.3207 + 0.0250i