Coding a function for the Fourier Transform.

I have the following code for a function here,

function F = FourierTransInv(vec,arg)
% This function caluclates the Discrete Fourier Transform or the Inverse
% Fourier Transform.
% The input vec is a vector of length 2^n (n in N with 0) where the entries 
% are all real numbers. arg is an input that determines whether or not the
% inverse or regular Fourier Transform is done. An arg of 1 corresponds to
% regular Fourier transform, while an arg of 0 is the inverse.
% The output F is the result of the Fourier or inverse transform.

if ~exist('vec')
    error('Please make sure your inputs are valid.');
elseif ~isrow(vec) || ~isnumeric(vec)
    error('Please make sure your inputs are valid.');
% Here 2^27 is the maximum because I observed that the difference in time
% between 2^20 and 2^21 is a factor of 2, where 2^20 took 30 seconds. So
% assuming the time increases by a factor of 2 each time, the minimum that 
% 2^27 can take is 1 hour which is an extremely long time.
% This can of course be changed easily.
elseif ~ismember(length(vec),2.^(0:27))
    error('Please make sure your input is a power of 2.');
elseif ~ismember(arg,[0,1])
    arg=1;
    warning('Argument is being set to 1 due to invalid input.');    
end

    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),arg);
            Fodd=FourierTransInv(vec(1:2:end),arg);
            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
        % Division by sqrt(N) because this part always divides twice.
        F=F/sqrt(N);
    end
end

Now this works for example if I input the vector and then run it back through the inversion I get the same vector back. However, I tried some more complicated things such as the one produced with this line of code,

n=2^8; L=6*pi; dt=L/n; t=dt*[0:n-1];
y = sin(5*2*pi*t/L)+sin(10*2*pi*t/L);

I've tested many more than this and it all seems to lead to the same conclusion, if I run it through my function, taking the Fourier Transform it agrees with the function fft( . ) everywhere which is great. Then once I run it back through I get back almost the same thing except that every term is multiplied by a factor of , why is this? I cannot see anywhere in the code where this would happen, even less so because the simple works but not some slightly more complicated vector.

Answer 1

[David Goodmanson]

Hi Anthony,

the only problem here is the normalization. For the inverse transform you are dividing by sqrt(N) every time, (N being the 'local' N). However, the ifft has an overall factor of 1/2^N where 2^N is the length of the initial vector that is being transformed. So you need to divide by 2 at each level of recursion, but not if you have an element of length 2^0 = 1. If you replace

if arg==0
    %  Division by sqrt(N) because this part always divides twice.
    F=F/sqrt(N);
end

with

if arg==0 & N~=1
    F=F/2;
end  

the results look good.