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

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?

• Why are you looping in Matlab? – mathreadler May 4 '16 at 17:07
• Apart from looping is very slow in matlab, the fourier transform can be normalized in different ways. This affects how the inverse transform will be. – mathreadler May 4 '16 at 17:19
• Compare your formula to the ones here en.wikipedia.org/wiki/Discrete_Fourier_transform that will be a good exercise in skimming and picking out relevant information too. – mathreadler May 4 '16 at 17:19
• So I guess the reason I am getting inconsistent answers is because of normalisation? – Ilya Lapan May 4 '16 at 18:41
• If everything looks right but there is a constant factor wrong on each sample I would say normalization is the most probable issue. Note that this is not necessarily "wrong", it just changes how you would need to implement the inverse transform. – mathreadler May 4 '16 at 18:54

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.