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)
采纳的回答
Abraham Boayue
2018-5-12
编辑：Abraham Boayue
2018-5-12
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
更多回答（3 个）
James Tursa
2018-4-2
编辑：James Tursa
2018-4-2
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 个评论
Julian Oviedo
2018-4-2
Erkan
2024-2-9
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
2018-5-11
1 个评论
Jan
2018-5-12
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.
Julian Oviedo
2018-5-11
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?
类别
在 帮助中心 和 File Exchange 中查找有关 Discrete Fourier and Cosine Transforms 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
