Can I make double discrete Fourier transform produce the reverse of data?

Question

Fourier[Fourier[{1,2,3,4}]] gives {1,4,3,2} instead of a mirror image {4,3,2,1}. Is there a way to modify the Fourier transform so that this additional shift does not happen?

The formula used for FFT is as follows:

{a, b} = {0, 1};
data = {1, 2, 3, 4};
n = Length[data]; s = Range[n];
f1 = 1/n^((1 - a)/2) Sum[ data[[r]] Exp[2 \[Pi]  I  b (r - 1) (s - 1)/n], {r, 1, n}]
f2 = 1/n^((1 - a)/2) Sum[ f1[[r]] Exp[2 \[Pi]  I  b (r - 1) (s - 1)/n], {r, 1, n}]

And it does produce the shift. Is it inevitable for discrete Fourier transform?

The question arose from my attempts to match the discrete Fourier transform of a sampled sequence with an analytical expression obtained for an infinite interval. None of the answers I tried gave a really accurate result. I did some guessing and came up with something like this:

dt = 0.2; L = 20;
pulse[t_] := Sin[t] Exp[-0.1 (t - 2)^2] + Exp[-t^2] - 0.5 Exp[-(t + 2)^2];
TS = Table[{t, pulse[t]}, {t, -L + dt, L, dt}];
fftshift[flist_] := RotateRight[flist, (Length@flist)/2 - 1];

FOURIER[LIST_] := Module[{dt, df, L, \[Mu], num, ord, timelist, ftlist},
dt = LIST[[2, 1]] - LIST[[1, 1]];
L = LIST[[-1, 1]] - LIST[[1, 1]];
num = Length[LIST];
\[Mu] = (L + dt)/Sqrt[2 \[Pi] (num)]; (* this +dt increases the accuracy significantly *)
df = 2 \[Pi]/(num  dt);
timelist = RotateLeft[LIST[[All, 2]], num/2 - 1];
ftlist = \[Mu] fftshift[Fourier[timelist, FourierParameters -> {0, 1}]];
Transpose[{df (Range[num] - num/2), ftlist}]
]
f[w_] = FullSimplify[FourierTransform[pulse[t], t, w]]
fft = FOURIER[TS];
Show[ListLinePlot[
Transpose[{fft[[All, 1]], Abs[fft[[All, 2]] - f[fft[[All, 1]]]]}], PlotRange -> {{-20, 20}, All}]]
ListLinePlot[{Transpose[{fft[[All, 1]], Re[fft[[All, 2]]]}], Transpose[{fft[[All, 1]], Re[f[fft[[All, 1]]]]}]}, PlotRange -> {{-10, 10}, All}]

This gives a very good accuracy of ~10^-15. However, when applied twice, this transformation doesn't produce f(-t), as it should in the case of the analytical Fourier, but f(-t-dt), where dt is the sampling step, just like the Fourier[] itself, so can it even be fixed?

ListLinePlot[{Transpose[{TS[[All, 1]], Reverse[TS[[All, 2]]]}], FOURIER[FOURIER[TS]]}, PlotRange -> All]

I have a feeling that something is missing.

Answer 1

Update:

I modified formula from Wikipedia link Generalized DFT.

sFT2[d_] := 
 With[{nn = Length[d], s = -(Length[d] - 1)/2}, 
  Total@Table[
    1/Sqrt[nn] d[[n]] Exp[-2 Pi I (k - 1 + s) (n - 1 + s)/nn], {n, 
     nn}, {k, nn}]]

isFT2[d_] := 
 With[{nn = Length[d], s = -(Length[d] - 1)/2}, (-1)^(nn + 1)
    Total@Table[
     1/Sqrt[nn] d[[n]] Exp[2 Pi I (k - s) (n - 1 + s)/nn], {n, 
      nn}, {k, nn}]]

Testing if it fulfills all requirements.

data = RandomInteger[{-10, 10}, RandomInteger[{1, 10}]];

(* testing if double fourier == reverse of data *)

Chop@N@Nest[sFT2, data, 2] == Reverse[data]

(* testing if inverse fourier of fourier == data *)

Chop@N@isFT2[sFT2[data]] == data

(* testing if double inverse fourier of double fourier == data *)

Chop@N@Nest[isFT2, Nest[sFT2, data, 2], 2] == data

---

True

True

True

Old version:

I took and rewrote definitions from this link:

Symmetry and periodicity of frequency-shifted discrete Fourier transform

sFT[d_, s_ : 0] := 
 With[{nn = Length[d]}, 
  Total /@ Table[
    1/Sqrt[nn] d[[n]] Exp[-2 Pi I (k - 1 + s) (n - 1)/nn], {k, 
     nn}, {n, nn}]]

isFT[d_, s_ : 0] := 
 With[{nn = Length[d]}, 
  Total /@ Table[
    1/Sqrt[nn] d[[n]] Exp[2 Pi I (k - 1) (n - 1 + s)/nn], {k, nn}, {n,
      nn}]]

sFT is discrete Fourier transform with frequency shift as second argument.

isFT is inverse discrete Fourier transform with frequency shift as second argument.

Here we test it with s=0.

data = Range@4;
s = 0;

sFT[data, s];
sFT[%, s]

isFT[%%, s] == data

---

{1, 4, 3, 2}

True

For s=0 it produces {1, 4, 3, 2}. (as usual discrete Fourier)

For s=1 it produces {4 I, -3, -2 I, 1}.

For s=2 it produces {3, -2, 1, -4}.

For s=3 it produces {2 I, 1, -4 I, -3}.

Unfortunately none of them is {4, 3, 2, 1}.

For non-integer s it produces complex output for double Fourier.

But maybe someone can modify the definitions so that the desired output of double Fourier is achieved.

Answer 2

The question is about the DFT (Discrete Fourier Transform) using Fourier[]. The Wolfram documentation under "Details and Options" states:

• The discrete Fourier transform $v_s$ of a list $u_r$ of length $n$ is by default defined to be $\frac1{\sqrt{n}}\sum_{r=1}^n u_r e^{2\pi i(r-1)(s-1)/n}.$
• Note that the zero frequency term appears at position 1 in the resulting list.

As the second bullet point explicitly states, the zero frequency term corresponds to position $1$. Also, the first bullet point formula has the factors $(r-1)(s-1)$ in the exponent corresponding to $u_r$ and $v_s$ which effectively means that the indexing should be regarded as zero based numbering. Also note that the indexing of the frequency transform is modulo $n$.

Using zero based numbering the definition is $ v_s = \frac1{\sqrt{n}}\sum_{r=0}^{n-1} u_r e^{2\pi i r s/n}. $ For the case $n=4$, note that $i = e^{2\pi i/4}$ and that the definition simplifies to $ v_s = \frac12\sum_{r=0}^3 u_r i^{r s}. $ Verify this with the code

u = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};
v = Fourier /@ u; w = Fourier /@ v;
Rationalize[2 v]
(* {{1, 1, 1, 1}, {1, I, -1, -I}, {1, -1, 1, -1}, {1, -I, -1, I}} *)
Rationalize[w]
(* {{1, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}, {0, 1, 0, 0}} *)

The result verifies that taking the Fourier transform twice gives $ w_s = \frac12\sum_{r=0}^3 v_r i^{r s} = u_{-s} = u_{4-s}. $

Note that $ w_0 = u_0 $ for all $n$, while $w_{n/2} = u_{n/2}$ if and only if $n$ is even.

Your original question was

Fourier[Fourier[{1,2,3,4}]] gives {1,4,3,2} instead of a mirror image {4,3,2,1}. Is there a way to modify the Fourier transform so that this additional shift does not happen?

If that is what you really want, then try this:

myFourier[v_] := With[{n=Length@v}, Fourier[v*I^(Range[n]/n*2)]*I^(Range[n]/n*2)*I^((n-3)/n)];
myIFourier[v_] := With[{n=Length@v}, InverseFourier[v/I^(Range[n]/n*2)]/I^(Range[n]/n*2)/I^((n-3)/n)];
Table[ Rationalize@Chop@Nest[myFourier, Range[n], 2], {n, 6}]
(* {{1}, {2, 1}, {3, 2, 1}, {4, 3, 2, 1}, {5, 4, 3, 2, 1}, {6, 5, 4, 3, 2, 1} *)

Which seems to do what you want as far as mirroring, but probably does not give useful Fourier transform.

The context for what you want to do is the dihedral group which Wikipedia states

In mathematics, a dihedral group is the group of symmetries of a regular polygon, which includes rotations and reflections.

In particular, you want a reflection symmetry which reverses the order of the elements of a list. The fixed point behavior of a reflection depends on the parity of the number $n$ of elements of the list. If $n$ is odd, then there is exactly one fixed point of every reflection. If $n$ is even, then there is either no fixed point or exactly two fixed points.

A previous result for zero frequency is $w_0 = u_0$ which is always a fixed point. This is not a problem if $n$ is odd, but if $n$ is even, then another fixed point is $w_{n/2} = u_{n/2}$. However, if $n$ is even, then a reversal of the list has no fixed points which is a contradiction. Thus, what you want is not possible for even $n$, but is possible for odd $n$ if the fixed point middle element is regarded as the zero frequency element.

Using the idea of shifting the data as explained in the Wikipedia DFT article, I wrote simple code which does what I think that you want for odd $n$ and is as close as possible for even $n$.

You can try using it as

sFT[u_] := With[{s = Quotient[Length@u, 2]}, RotateRight[Fourier@RotateLeft[u, s], s]];
isFT[u_] := With[{s = Quotient[Length@u, 2]}, RotateRight[InverseFourier@RotateLeft[u, s], s]];
Table[ Rationalize@Chop@Nest[sFT, Range[n], 2], {n, 6}]
(* {{1}, {1, 2}, {3, 2, 1}, {1, 4, 3, 2}, {5, 4, 3, 2, 1}, {1, 6, 5, 4, 3, 2}} *)