1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Why are you looping in Matlab? $\endgroup$ – mathreadler May 4 '16 at 17:07
  • $\begingroup$ 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. $\endgroup$ – mathreadler May 4 '16 at 17:19
  • 1
    $\begingroup$ 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. $\endgroup$ – mathreadler May 4 '16 at 17:19
  • $\begingroup$ So I guess the reason I am getting inconsistent answers is because of normalisation? $\endgroup$ – Ilya Lapan May 4 '16 at 18:41
  • $\begingroup$ 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. $\endgroup$ – mathreadler May 4 '16 at 18:54
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.