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.
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.
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.
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
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 ...
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
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--');
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
"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
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!
