3
$\begingroup$

Edit 1/30 - Taking @Fat32's edits into account, it seems like there is still an issue with the scale of the frequency axis. While version 2 correctly identifies the response at 1 HZ, version 1 seems to be covering a frequency range of 1/T.

I'd like version 1 to cover the same range of frequencies as version 2. Can anyone help with this?

enter image description here

(End edit 1/30)

--

I have 2 DFT algorithms that I would expect to return same frequency response. Unfortunately, the more efficient recursive algorithm does not seem to behave as expected.

Version 1 - Efficient Divide and Conquer method, but output is incorrect:

function [ s ] = DFT_ver_1( x )

N = length(x);

if N == 1
    s = x(1);
    return;
else

X1 = DFT_ver_1(x(1:2:N));
X2 = DFT_ver_1(x(2:2:N)); 
s = zeros(N,1);

for k = 1:(N/2)
        W(k)     = exp(-1j*2*pi*(k + 1)/N);
        s(k)     = X1(k) + W(k) * X2(k);
        s(k+N/2) = X1(k) - W(k) * X2(k);  
end

Version 2 - Inefficient, but generates expected results:

function [ s ] = DFT_ver_2( x , f, n)

N = length(x);

b = zeros(N,N);
b_1 = -2*pi*(f/N);

for idx = 1 : N
for jdx = 1 : N
    b(jdx, idx) = b_1(idx) * n(jdx);
end
end

anlyz_fn = cos(b) + 1i * sin(b);
s = single(zeros(N,1));

for k = 1 : N
    val = 2 * sum(anlyz_fn(k,:) * x.');
    s(k) = val;
end

end

I would expect both algorithms to return the same result. However, this does not appear to be the case.

You can see that the bottom plot (version 2) shows the correct frequency response to a 1 HZ signal, but the middle (version 1) does not:

enter image description here

Code to run both:

       T  = single(3.20); % (sec) time window
       dt = single(0.05); % (sec) sample time
        N = T / dt;          % (ct) num samples
        n = 0 : N - 1;       % (ct) bucket index vector
        t = 0 : dt : T - dt; % (sec) time vector
smpl_freq = 1 / dt;          % (Hz) sampling frequency
 freq_res = smpl_freq / N;   % (Hz) frequency resolution
        f = n * freq_res;    % (Hz) freq. bucket vector

x_t = sin( 2*pi*t );

subplot(3,1,1);
plot(t, x_t);
xlabel('Sample Time (sec)');
ylabel('Magnitude');

subplot(3,1,2);

s_v1     = DFT_ver_1( x_t); %, f, n );
s_v1_mag = abs( s_v1(1:N/2) )/N; 

stem(f(1:N/2)/T, s_v1_mag);

xlabel('Freq. Response (Hz) v1');
ylabel('Magnitude');

s_v2     = DFT_ver_2( x_t , f, n );
s_v2_mag = abs( s_v2(1:N/2) )/N; 

subplot(3,1,3);
stem(f(1:N/2) / T , s_v2_mag);
xlabel('Freq. Response (Hz) v2');
ylabel('Magnitude');
$\endgroup$
  • 1
    $\begingroup$ Hey, just a remark: In principle, the recursive one seems to be an implementation of the Cooley-Tukey algorithm of the FFT. On a machine that can run recursion as fast as iterative algorithms, it would use fewer computations for the same result. In practice a) Matlab is really inefficient at handling its own programming language, so the fact you're doing this with a for loop instead of with a stupid DFT matrix multiplication means this is way way slower than the "obvious" linear algebra method of considering the DFT a change of basis, for anything but humongous DFTs, $\endgroup$ – Marcus Müller Jan 29 at 22:25
  • $\begingroup$ and b) Matlab is not optimized for recursion, so don't expect either version to be reliably faster or slower than the other. Both are slow. Nobody would realistically implement an FFT in Matlab itself. Matlab's fft function uses an external library, FFTw. $\endgroup$ – Marcus Müller Jan 29 at 22:26
  • 1
    $\begingroup$ @MarcusMüller Nobody would realistically implement an FFT in Matlab itself. ... well I've implemented (decades ago) my first FFT routine in MATLAB :-)) Isn't it the first thing to do in DSP: implement it in MATLAB ;-)) $\endgroup$ – Fat32 Jan 29 at 23:28
  • 2
    $\begingroup$ @Fat32 Right, but you do it once, then use the built-in commands :) Actually, I like to use my own DTFT implementation in Matlab/Julia. For the applications I do, the transform is seldom the bottleneck. $\endgroup$ – MBaz Jan 29 at 23:53
  • 1
    $\begingroup$ @Fat32 we've all been there :) Yeah, I did an "FFT" in matlab, too. We just this year scrapped that exercise from our digital comms computer lab; it was the "ok, you've heard of the FFT. Here's C-T and stubs, go forth and implement it" exercise to get the students to warm up to Matlab coding and get them to see the n·log n complexity. Coding took ages, computation took ages, and the fact that "compute the DFT matrix and use it to compute many transforms" was much faster than student's "F"FT proved that we were trying to show too many things at once. $\endgroup$ – Marcus Müller Jan 30 at 8:21
1
$\begingroup$

All right, you made two small mistakes and here is the corrected version: (I changed your function's name). Btw, you can get rid of the for loop too... ;-)

function [ s ] = recDFT2( x )
%
% a recursive DFT
%


N = length(x);

if N == 1
    s = x(1);
    return;
else

    X1 = recDFT2(x(1:2:N));
    X2 = recDFT2(x(2:2:N)); 
    s = zeros(N,1);

    for k = 1:(N/2)
        W(k)     = exp(-1j*2*pi*(k - 1)/N);
        s(k)     = X1(k) + W(k) .* X2(k);
        s(k+N/2) = X1(k) - W(k) .* X2(k);  
    end
end
$\endgroup$
  • $\begingroup$ ty @Fat32, my results are starting to make more sense. It seems now that the frequency axis range of version 1 is scaled by a factor of 1/T to version 2. Could you check out the edit that I made? $\endgroup$ – darkpbj Jan 30 at 21:35
  • $\begingroup$ the so called Frequency Axis Scaling is not related with DFT. It just works fine now. Gets in a sequence of x[n] and produces the correct output X[k]. That scale is about how you display the result, other part of the code. And it seems you make another mistake there but that's not about DFT... $\endgroup$ – Fat32 Jan 30 at 21:39
  • $\begingroup$ Thanks, I can see how that is the case, but could you explain why the two algorithms are producing results at two different frequency resolutions? $\endgroup$ – darkpbj Jan 31 at 19:37

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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