FFT implementation by myselft
10 views (last 30 days)
Show older comments
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)
0 Comments
Accepted Answer
Abraham Boayue
on 12 May 2018
Edited: Abraham Boayue
on 12 May 2018
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
More Answers (3)
James Tursa
on 2 Apr 2018
Edited: James Tursa
on 2 Apr 2018
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).
2 Comments
Erkan
on 9 Feb 2024
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?
Julian Oviedo
on 11 May 2018
1 Comment
Jan
on 12 May 2018
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.
See Also
Categories
Find more on Signal Processing Toolbox in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!