2
$\begingroup$

I am implementing DFT in Octave. Here's my code:

function [real_comp, imag_comp] = mydft (samples)
    N = columns(samples);
    for m = 1:N
        real_accum = 0;
        imag_accum = 0;
        for n = 1:N
            real_factor = cos((2 * pi * (n-1) * (m-1)) / N);
            # imag_factor = -I * sin((2 * pi * (n-1) * (m-1)) / N);
            imag_factor = sin((2 * pi * (n-1) * (m-1)) / N);
            real_accum += real_factor * samples(n);
            imag_accum += imag_factor * samples(n);
        endfor
        real_comp(m) = real_acuum;
        imag_comp(m) = imag_accum;
    endfor
endfunction

My goal is to compute the magnitude at each DFT bin and plot it. I am passing in a list of real-only samples and storing the real and imaginary components in different variables. Here's the sample vector that I am using:

samples = [0.3535, 0.3535, 0.6464, 1.0607, 0.3535, -1.0607, -1.3535, -0.3535]

The samples are carefully chosen to prevent spectral leakage at fs = 8KHz. The graph of the DFT magnitude is exactly what I expect. There's nothing wrong with that. But what I don't get is that some components in my 8 point DFT which are supposed be zero, are not exactly zero but rather a very small, close-to-zero value. The built-in fft() function, on the other hand, has zeroed out components where I expect.

For example, I print the elements of output vector produced by my DFT and octave's FFt using the code below:

fft_out = fft(samples);
[real, imag] = mydft(samples);
dft_out = real + (-I * imag);
for i=1:columns(samples)
    disp(fft_out(i));
    disp(dft_out(i));
    printf("\n");
endfor

The point-to-point comparison is given in the table below:

bin_index fft_out dft_out
1 -1.0000e-04 -1.0000e-04
2 0.00000 - 3.99988i 4.7184e-16 - 3.9999e+00i
3 1.4141 + 1.4144i 1.4141 + 1.4144i
4 0.0000e+00 - 8.0820e-05i 1.3878e-16 - 8.0820e-05i
5 -1.0000e-04 -1.0000e-04 - 1.4350e-16i
6 0.0000e+00 + 8.0820e-05i 2.3037e-15 + 8.0820e-05i
7 1.4141 - 1.4144i 1.4141 - 1.4144i
8 0.00000 + 3.99988i -2.0817e-15 + 3.9999e+00i

Now, I wouldn't care too much about the tiny difference but if I compute the phase angles of fft_out and dft_out,

disp(arg(fft_out))
disp(arg(dft_out))

I get the sign flipped in some cases:

fft_out: 3.14159  -1.57080   0.78550  -1.57080   3.14159   1.57080  -0.78550   1.57080
dft_out:-3.14159  -1.57080   0.78550  -1.57080  -3.14159   1.57080  -0.78550   1.57080

So, is there something wrong with my DFT implementation? Is it just a rounding or precision issue?

$\endgroup$
1
  • 1
    $\begingroup$ My my. That 'samples' input sequence certainly does look familiar! $\endgroup$ – Richard Lyons Jan 3 at 11:24
3
$\begingroup$

These are just rounding errors, and you only get a flipped sign for phase values of $\pm\pi$, which is not an issue because the phase value is ambiguous, i.e., adding or subtracting multiples of $2\pi$ doesn't change the value of the complex coefficient.

I assume that your code is just for an educational purpose, because it is neither efficient nor numerically accurate.

$\endgroup$
2
  • $\begingroup$ Yes. It is just an exercise. But I would still like to know what numerical imperfections do you see :) $\endgroup$ – Devashish Jaiswal Jan 3 at 11:16
  • 2
    $\begingroup$ @DevashishJaiswal: The problem is the accumulation of numbers, which can become inaccurate if the magnitudes of those number vary strongly. The inaccuracies come mainly from the fact that you need many more operations than an FFT algorithm would need. Also, Matlab is vector/matrix-based, so nested for-loops are highly inefficient and should (and can) be avoided. $\endgroup$ – Matt L. Jan 3 at 11:51
1
$\begingroup$

Since you didn’t (and can’t, unless you use strictly symbolic transformations) use infinite precision arithmetic and values, you are seeing numerical noise instead of zeros (etc.). And the phase of noise values can be random.

An FFT does less arithmetic than a DFT, so has less opportunities to add more numerical noise.

$\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.