FFT implementation by myselft

14 次查看(过去 30 天)
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
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
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 个评论
Erkan
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
Julian Oviedo 2018-5-11
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.
  1 个评论
Jan
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
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?

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by