Discrete Fourier transform implementation giving results that are order of magnitude off

Question

I tried implementing a Discrete Fourier transform in Matlab, but I found my results an order of magnitude off. I used next definition of DFT: $$ F(u) = \frac{1}{2N} \sum^{N-1}_{x=-N} f(x) e^{- \pi i x u/N} $$ I obtained separately the values for real and imaginary coefficients from this expressions and implemented the algorithm in Matlab. I set N=500, so I made an array of x values of size 500, and array of the result values for real and imaginary of size 500. Here is my Matlab code of my DFT function:

function [ rFT, iFT ] = DFT( N, rf, imf )

rFT = zeros(length(rf),1);
iFT = zeros(length(imf),1);

a = fix(length(rf)/2) + 1;
for u=(-N):(N-1) 
    sumR = 0;
    sumI = 0;
    for x=(-N):(N-1) 
        sumR = sumR + rf(x+a) * cos(- pi*(x)*u/N) - imf(x+a)*sin(- pi*(x)*u/N); %adjust for the fact that something starts counting at 0, not -250 
        sumI = sumI + rf(x+a) * sin(- pi*(x)*u/N) - imf(x+a)*cos(- pi*(x)*u/N);

    end 
    rFT(u+a) = sumR;
    iFT(u+a) = sumI;
end
rFT = rFT/(2*N);
iFT = iFT/(2*N);
end

Now, if I try to input a box function for example of width 10 and height 1 centered on origin, I find that the result looks like a sinc function, which is what I expected but the maximum value of the function is of order $10^{-2}$ (value at $u=0$), while I find just by taking the integral myself that the max should be 10. And similarly, if I take a transform of a convolution of the box with itself (I calculated the convolution on paper) and compare it to the square of the transform of the box, I find they look similar in shape, but again off by a very big scalar constant. Is my understanding of theory wrong, or is there something off with my implementation?

Answer 1

It seems to me that your implementation is correct but that you use a non-unitary transform: that is you need different normalization factors for forward and backward transform. You can for example turn your DFT into a unitary DFT by replacing

rFT = rFT/(2*N);
iFT = iFT/(2*N);

with

rFT = rFT/(2*N).^0.5;
iFT = iFT/(2*N).^0.5;

However in the Matlab programming language using loops is very slow compared to doing matrix and vector operations. By comparison running DFT on a 127 sample long $\sin^2(x)$ function ( in seconds, and actually I run it in Octave, but it is designed to be compatible ):

t_vec =  0.0024350
t_for =  1.1651
t_for/t_vec
ans =  478.49

A factor of almost 500 times! So if you expect to use matlab or other vectorized languages in the future for more involving and computationally intensive things I would encourage you to learn thinking and writing code which takes advantage of that.