I have re-created a Matlab code to compute DTFT of a signal. The example is from a text book.
The following is the code:
n = 1:10;
x = (0.9*exp(j*pi/3)).^n;
k = -200:200;
w = (pi/100)*k;
b_var = (exp(-j*pi/100)).^(n'*k);
X = x*b_var; magX = abs (X);
angX = angle(X);
realX = real(X);
imagX = imag (X);subplot (2,1,1);
plot (w/pi, magX); grid
xlabel ('frequency in pi units'); ylabel('|X|')
title ('Magnitude part');subplot (2,1,2);
plot (w/pi, angX/pi); grid
xlabel ('frequency in pi units'); ylabel('radians/pi')
title ('phase part');
The example works as shown in the text book. However, I am unable to understand the following line:
b_var = (exp(-j*pi/100)).^(n'*k);
Why exponent raised by (n'*k) is being computed.
As far I know, DFT is computed as X(k) = summation(exp(-j*2*pi*k*n/N))
So, exp((-j*2*pi*k*n)/N)) is expected but I am confused what exp(-j*pi/100)^(n'*k) does.
Note: I have ignored "." operator in ".^" just to simply the understanding.
![](https://d23s79tivgl8me.cloudfront.net/user/profilepictures/14446.jpg)
Hi Sharan123. Keep in mind, the DTFT is an equation and the 'w' (omega) frequency variable in that equation is continuous. We cannot use or "compute" continuous variables in a computer using software. What we CAN do is evaluate the DTFT at various "fixed" values of 'w', and doing that is called the discrete Fourier transform (DFT).
I wonder, what does the phrase "re-created" mean in your first sentence? Did you copy someone's code or did your modify someone's code? Why does your frequency variable 'w' range from w = -2pi -to- +2pi? That seems a bit strange to me.
I think that code is equivalent to computing a 401-point DFT (over a 4pi freq range) of your damped 'x' complex sinusoid that has been zero-padded out to 401 samples. And I think the code's author is being super clever in using MATLAB's matrix capabilities by employing that "b_var = (exp(-j*pi/100)).^(n'*k);" command which, somehow, eliminates the need to actually zero pad your 'x' input. But I'm too lazy to experiment to figure out exactly what that "b_var = (exp(-j*pi/100)).^(n'*k);" command is doing. What fascinates me is that I do not see a summation operation in your code, but the code DOES seem to being computing a correct "zero padded" DFT!
Sorry that I could be of such little help.
![](https://www.embeddedrelated.com/new/images/defaultavatar.jpg)
Dear Rick,
You are right about computing DTFT at 401 discrete points.
The idea is to evaluate X(e^j*omega) and examine its periodicity. That's the reason why it is being computed from -2*pi to 2*pi
The lack is summation puzzled me too but I have analyzed it step by step and things look fine. I will break it down further and post more (on MOnday) explaining what is happens at each step. But right now I am stuck at b_var.
Anyway, thanks a lot. I will spend some more time and post my understanding back on this forum.
![](https://d23s79tivgl8me.cloudfront.net/user/profilepictures/37480.jpg)
The line: x = (0.9*exp(j*pi/3)).^n;
generates a distorted single tone test vector
The line: b_var = (exp(-j*pi/100)).^(n'*k);
generates base tones across frequency grid (10 points x 401 samples)
The line: X = x*b_var;
correlates input with base functions using matrix multiplication
This correlation is faster than vector summation.
The frequency axis is not quite right being -2pi~+2pi
![](https://www.embeddedrelated.com/new/images/defaultavatar.jpg)
Kaz,
This is DFT computation. I am not able to understand where the term (exp(-j*2*pi*k*n/N)) is computed in the DFT expression is computed. Obviously, something of the form (a+jb)^k is being done but how exp(-j*2*pi*k*n/N) is decomposed into this, I am not able to understand ...
![](https://d23s79tivgl8me.cloudfront.net/user/profilepictures/37480.jpg)
It is DFT.
DFT is based on correlating input with cos/sin base functions that are generated over the chosen frequency grid.
There are three steps in this exercise:
1) generate an input
2) generate base functions (as set of cos/sin across grid)
3) correlate base functions with input (summation or matrix)
I am no fan of matrix computations but it does make it quicker.
Below is my code using loops to do DFT (very slow).
Then I compare it with matlab fft and the speed is quite different.
I am pretty sure Matlab uses matrix instead of loop.
--------------------------------------------------
n = 2^8;
x = randn(1,n);
y = zeros(1,n);
for k = 1:n
for i = 1:n
y(k) = y(k) + x(i)*exp(-j*2*pi*(i-1)*(k-1)/n);
end
end
y = round(y*32767);
%compare to fft function, both rounded
ref = fft(x,n);
ref = round(ref*32767);
plot(y - ref);
---------------------------------------------------
Below is alternative direct way of generating base functions for each and all values of f in one single step:
n = 31;
f = linspace(-.5,.5,n)'; %frequency grid
F = exp(-j*2*pi*f*[-(n-1)/2:(n-1)/2]); %grid of sinusoids
![](https://d23s79tivgl8me.cloudfront.net/user/profilepictures/37480.jpg)
My guess is your code is generating single tone (x) & base functions wrong and your output (X) does not compare well with matlab fft. This apart from wrong frequency axis.
Here is direct comparison between matlab fft & my matrix approach:
n = 1005;
x = exp(j*2*pi*(0:n-1)*.01); %test tone
%matlab fft
test = fft(x);
%matrix:
f = linspace(0,1,n)'; %frequency grid
F = exp(-j*2*pi*f*[0:n-1]); %grid of cos/sin
X = x*F;
plot(abs(test));
hold
plot(abs(X),'r--');
![](https://d23s79tivgl8me.cloudfront.net/user/profilepictures/103185.jpg)
Sharan123,
The DFT *is* a matrix multiplication.
Z = F * S
Where Z is your DFT bin set, F is a matrix composed of the sinusoidal basis vectors, and S is your signal.
b_var = (exp(-j*pi/100)).^(n'*k);
Constructs a truncated version of F. Since your signal only has 10 samples, there is no point in using the full array. This is the zero padded equivalency.
X = x*b_var;
Does the DFT.
See:
https://en.wikipedia.org/wiki/DFT_matrix
Hope this helps,
Ced
![](https://d23s79tivgl8me.cloudfront.net/user/profilepictures/103185.jpg)
"something of the form (a+jb)^k is being done but how exp(-j*2*pi*k*n/N) is decomposed into this"
exp(-j*2*pi*k*n/N) = e^(-j*2*pi*k*n/N)
= [e^(-j*2*pi*n/N)]^k
= [cos(2*pi*n/N)-j*sin(2*pi*n/N)]^k
a = cos(2*pi*n/N)
b = -sin(2*pi*n/N)
Look up Euler's equation, or check out my first blog article.
Hope this helps, too.
Ced
![](https://d23s79tivgl8me.cloudfront.net/user/profilepictures/14446.jpg)
Hi Sharan123. kaz's May 13 reply made me realize my brain was not working when I posted my above reply. No explicit summation command in software is needed because of the matrix multiplications. The product of one row of a matrix by a column matrix implicitly results in a summation of individual products. The summation is "built-in" during the row-by-column multiplication!