radix 2 disspation in frequency

OK

[Walter Roberson]

That image is too blurry for me to read. You should post your code instead of an image of the code.

[ZHAO HONG QIU]

This is my file. Thanks.

Answer 1

[Geoff Hayes]

Zhao - I think the problem may be due to trying to calculate the even and odd terms in your inner for loop (perhaps to improve computation time). If I remember correctly, the DFT (without a shift) can be simplified to just

    X = zeros(N,1);
    for k=0:N-1
        innerSum = 0;
        for n=0:N-1
            innerSum = innerSum + x(n+1)*(cos(2*pi*n*k/N) - 1i*sin(2*pi*n*k/N));
        end
        X(k+1) = innerSum;
    end

where X is your sum_all. But note how as we iterate from 0 to N - 1 in the inner for loop, that the factor x(n+1) is used but use n within the cos and sin function calls. That is not the case in your code where we have factors of 2n or 2n-1 but use just n within the cos and sin. I also don't understand why divide the N by two or the other factor for the sum_odd(k) terms of

 (cos(2*pi*(k-1)/N)-1j*sin(2*pi*(k-1)/N))

Please clarify why these extra steps are being made.

If you do use the above simplified DFT code, then we will observe the frequency at 300 Hz if we do

 f = (0:N-1)*(Fs/N);
 subplot(2,1,2);
 plot(f,abs(X)/(N/2));

which is still very similar to what you have written.

[ZHAO HONG QIU]

HI Geoff, I want to achieve radix-2 FFT algorithm. Reference http://www.cmlab.csie.ntu.edu.tw/cml/dsp/training/coding/transform/fft.html

[Geoff Hayes]

Hi Zhao - I see what you are trying to do with the above. I think the problem is that there is some confusion over what are the even terms and which are the odd ones. This is especially difficult when translating from an algorithm which assumes zero-based arrays into MATLAB which uses one-based arrays. In the above equations, the even terms would correspond to

 x[0], x[2], x[4], ..., x[N-2]

with the odd terms as

 x[1], x[3], x[5], ..., x[N-1]

In MATLAB, the even terms would then be

 x(1), x(3), x(5), ..., x(N-1)

with the odd terms as

 x(2), x(4), x(6), ..., x(N)

Your code would then be (if we continue with iterating from 1 to N/2)

for k=1:N/2
   Xre=0; %each one always reset the sum 
   Xim=0;
   Xre1=0; 
   Xim1=0;
   for n=1:N/2
       %DFT even of x[n]
       Xre = Xre+x(2*n)*cos(2*pi*(k-1)*2*(n-1)/N);
       Xim = Xim+x(2*n)*(-1j)*sin(2*pi*(k-1)*2*(n-1)/N);
       %DFT odd of x[n]
       Xre1 = Xre1+x(2*n-1)*cos(2*pi*(k-1)*2*(n-1)/N);
       Xim1 = Xim1+x(2*n-1)*(-1j)*sin(2*pi*(k-1)*2*(n-1)/N);
   end 
   sum_even(k) = (cos(2*pi*(k-1)/N)-1j*sin(2*pi*(k-1)/N))*(Xre+Xim); 
   sum_odd(k)  = (Xre1+Xim1);
   sum_all(k) = sum_even(k) + sum_odd(k);
end

Note how I have replaced the (n-1)/(N/2) with 2(n-1)/N so that the code more closely resembles the equation. I also removed the final division by N because it doesn't appear in the above algorithm and is left for when we plot the data as

 plot(f,abs(sum_all)/(N/2));

When coding an algorithm, it is sometimes better to stay true to the algorithm (set of equations) that you are trying to implement so that it is easier to avoid making mistakes or it may make it clear where a mistake or problem lies. In this case, the above equations could simplify to

sum_even = zeros(N/2,1);
sum_odd  = zeros(N/2,1);
sum_all  = zeros(N/2,1);
for k = 0:N/2-1
   for n=0:N/2-1
       %DFT even of x[n]
       sum_even(k+1) = sum_even(k+1) + x(2*n+1)*(cos(2*pi*k*2*n/N) ...
                                     +     (-1j)*sin(2*pi*k*2*n/N));
       %DFT odd of x[n]
       sum_odd(k+1) =  sum_odd(k+1)  + x(2*n+2)*(cos(2*pi*k*(2*n)/N) ...
                                     +     (-1j)*sin(2*pi*k*(2*n)/N));
   end 
   sum_all(k+1) = sum_even(k+1) + sum_odd(k+1)*(cos(2*pi*k/N) - 1j*sin(2*pi*k/N));
end
f = (0:N/2-1)*(Fs/N);
subplot(2,1,2)
plot(f,abs(sum_all)/(N/2));

[ZHAO HONG QIU]

Thank you very much! It finally was successful!