Cooley-Tuckey Algorithm, Inverse and Forward.

Question

I'm trying to make some code in MatLab that takes a vector of length $2^n$ (a sequence) and spits out it's Discrete Fourier Transform, or it's inverse Discrete Fourier Transform depending on another input. Now I have the regular Discrete Fourier Transform working, however I am stuck on the inverse. After messing around a bit I noticed that it seems like the Cooley-Tuckey Algorithm seems to work in reverse with some minor differences. So suppose I have a sequence $\{F_k\}_{k \in \mathbb{N}}$, and I want to find it's inverse $f_k$, then I'd use the following,

$$f_k=\frac{1}{N}\sum_{j=1}^{N}F_j\exp\left(\frac{2\pi i}{N}(j-1)(k-1)\right).$$

The $k-1$ because I want to start from $f_1$ as opposed to $f_0$. Which is essentially the same as the regular Discrete Fourier Transform without the minus sign and an extra $1/N$ factor. As such I attempt to do the following ( with $\omega_N=\exp(2\pi i/N)$ ).

$$I_k :=Nf_k=\sum_{j=1,\text{ odd}}^{N}F_j\omega_N^{(j-1)(k-1)}+\sum_{j=1,\text{ even}}^{N}F_j\omega_N^{(j-1)(k-1)}.$$

$$I_k=\sum_{j=1}^{N/2}F_{2j-1}\omega_N^{(2j-2)(k-1)}+\sum_{j=1}^{N/2}F_{2j}\omega_N^{(2j-1)(k-1)}.$$

$$I_k=\sum_{j=1}^{N/2}F_{2j-1}(\omega_N^2)^{(j-1)(k-1)}+\omega_N^{k-1}\sum_{j=1}^{N/2}F_{2j}(\omega_N^2)^{(j-1)(k-1)}.$$

$$I_k=\sum_{j=1}^{N/2}F_{2j-1}(\omega_{N/2})^{(j-1)(k-1)}+\omega_N^{k-1}\sum_{j=1}^{N/2}F_{2j}(\omega_{N/2})^{(j-1)(k-1)}.$$

Then we are basically done because we have $I_k$ in terms of Fourier Transforms of the odd and even sequences. Therefore, the inverse is just the same as the Cooley-Tuckey Algorithm except with $\omega_N$ and dividing every entry by $1/N$ should get you the result you want.

Then if this is correct then the only problem is in my code, though I'm not entirely sure if I should be asking such a question here, this is the code I'm using.

function F = FourierTransInv(vec,arg)
     N=length(vec);
     F=ones(1,N);   
     switch arg
         case 1
             w=exp(-2*pi*1i/N);
         case 0
             w=exp(2*pi*1i/N);
     end
     switch N
         case 1         
             F(1)=vec(1);
         otherwise
             Feven=FourierTransInv(vec(2:2:end),1);
             Fodd=FourierTransInv(vec(1:2:end),1);
             for k=1:N/2
                F(k)=Fodd(k)+w^(k-1)*Feven(k);
                F(k+N/2)=Fodd(k)-w^(k-1)*Feven(k);
             end
     end    
     if arg==0
        F=F/N;
     end
end

As such I've implemented that if we want the Inverse then simply $\omega_N$ is as above and at the end $F$ is divided by $N$. The strange thing is when I run this code, every odd entry agrees with the true result we want, e.g when I run the Fourier Transform of $[1,2,3,4,5,6,7,8]$, then run it back through the inverse, it agrees with $1,3,5,7$ only, and the rest are just gibberish.

Just for reference the output for this input specifically is, $[1+0i,6+2i,3+0i,4-2i,5+0i,6-2i,7-0i,4+2i]$.

If the regular Fourier Transform is working as intended (I tested it for many sequences using online calculators and they all seem to agree), I don't see a reason why the inverse shouldn't work assuming what I did before showing the code is correct, since I basically reuse what was already working.

Answer 1

I don't know whether there are other issues, but here is one bug that I find.

In the recurrence calls Feven = ... and Fodd = ..., you should pass arg instead of 1.

That is,

Feven=FourierTransInv(vec(2:2:end),arg);
Fodd=FourierTransInv(vec(1:2:end),arg);