FFT implementation by myselft

Hello. Im studying discrete fourier transform and for that I want to implement it. I mean, this is the code I made but there is a problem

I defined a little discrete signal, x. And I calculate the DFT and then the inverse DFT but I dont get the same signal.

What Im doing wrong?

Thanks in advance.

    x=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1] 

sum=0
for k=1:40 

     for j=1:40

       sum= sum+x(j)*exp(- 1i * 2 * pi * (1/20)*(j-1)*(k-1))  

     end

     a(k)=sum
      sum=0
  end

n=1:1:40 

sum=0
for k=1:40 

     for j=1:40

       sum= sum+a(j)*exp( 1i * 2 * pi * (1/20)*(j)*(k))  

     end

     b(k)=(1/40)*sum
      sum=0
  end

subplot(1,2,1)
plot(n, x)
subplot(1,2,2)
plot(n,b)

Answer 1

[Abraham Boayue]

It would be more efficient to use a vector form of the DFT for the coding in matlab. As an example, see below. Note that this method is only best for shorter sequences ,x(n), if you have a larger value for N, then another method would be preferred.

clear variables
close all
%%Implementing the DFT in matrix notation
xn = [0 2 4 6 8];        % sequence of x(n)
N = length(xn);        % The fundamental period of the DFT
n = 0:1:N-1;           % row vector for n
k = 0:1:N-1;           % row vecor for k
WN = exp(-1i*2*pi/N);  % Wn factor 
nk = n'*k;             % creates a N by N matrix of nk values
WNnk = WN .^ nk;       % DFT matrix
Xk = xn * WNnk;        % row vector for DFT coefficients
disp('The DFT of x(n) is Xk = ');
disp(Xk)
magXk = abs(Xk);       % The magnitude of the DFT
%%Implementing the inverse DFT (IDFT) in matrix notation
n = 0:1:N-1;           
k = 0:1:N-1;           
WN = exp(-1i*2*pi/N);  
nk = n'*k;             
WNnk = WN .^ (-nk);    % IDFS matrix
x_hat = (Xk * WNnk)/N; % row vector for IDFS values
disp('and the IDFT of x(n) is x_hat(n) = ');
disp(real(x_hat))
% The input sequence, x(n)
subplot(3,1,1); 
stem(k,xn,'linewidth',2); 
xlabel('n'); 
ylabel('x(n)')
title('plot of the sequence x(n)')
grid
% The DFT
subplot(3,1,2); 
stem(k,magXk,'linewidth',2,'color','r'); 
xlabel('k'); 
ylabel('Xk(k)')
title('DFT, Xk')
grid
% The IDFT
subplot(3,1,3); 
stem(k,real(x_hat),'linewidth',2,'color','m'); 
xlabel('n'); 
ylabel('x(n)')
title('Plot of the inverse DFT, x_{hat}(n) = x(n)')
grid

[Vaibhav Pawale]

thanks sir!!!

Answer 2

[James Tursa]

Use 1/40 instead of 1/20, and modify your inverse formula to account for the -1 difference between the MATLAB indexing (1-based) and the indexing used by the inverse formula (0-based). E.g.,

     for j=1:40
       sum= sum+x(j)*exp(- 1i * 2 * pi * (1/40)*(j-1)*(k-1))  
     end

and

     for j=1:40
       sum= sum+a(j)*exp( 1i * 2 * pi * (1/40)*(j-1)*(k-1))  
     end

It would be better to use a variable name other than "sum", since that is the name of an intrinsic MATLAB function. Also, it would be better to generalize your code using variables for indexing limits instead of hard-coded numbers (e.g., N instead of 40).

[Julian Oviedo]

tanks a lot. It did work

[Erkan]

Hi James, how can we calculate the w(2*pi*f) in the jw term in the Fourier integral formula from this

(- 1i * 2 * pi * (1/40)*(j-1)*(k-1)) expression? So how can we find frequency values?

Answer 3

[Julian Oviedo]

Hello. Ive other question regards this topic.

if the input sinal is not periodic. I have to take out the N=40?

Can it be:

if true
 sum=0
for k=1:40 
   for j=1:40
     sum= sum+x(j)*exp(- 1i * 2 * pi *(j-1)*(k-1))  
   end
   a(k)=sum
    sum=0
end
n=1:1:40 
sum=0
for k=1:40 
   for j=1:40
     sum= sum+a(j)*exp( 1i * 2 * pi *(j)*(k))  
   end
   b(k)=(1/40)*sum
    sum=0
end
subplot(1,2,1)
plot(n, x)
subplot(1,2,2)
plot(n,b)
end

Thanks in advance.

[Jan]

Please ask one question per thread only. Injecting a new question in the section for answers is rather confusing and it is not longer clear, to which question the answers belong.

If James' answer solves your problem, please accept it. Open a new thread for further questions. Thanks.

Answer 4

[Julian Oviedo]

Because I have an array of 1024 values and I make:

   sum=0
    for k=1:1024 
     for j=1:1024
       sum= sum+x(j)*exp(- 1i * 2 * pi *(j-1)*(k-1))  
     end
     a(k)=sum
      sum=0
    end

and what I got (amplitude not phase) is this (the frecuency are bad!)

What Im doing wrong?