Forums

DTFT of a signal

Started by Sharan123 3 years ago9 replieslatest reply 3 years ago535 views

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.

[ - ]
Reply by Rick LyonsMay 12, 2017

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.




  

[ - ]
Reply by Sharan123May 13, 2017

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.

[ - ]
Reply by kazMay 13, 2017

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 

[ - ]
Reply by Sharan123May 13, 2017

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 ...

[ - ]
Reply by kazMay 13, 2017

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

[ - ]
Reply by kazMay 13, 2017

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--');

[ - ]
Reply by CedronMay 13, 2017

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

[ - ]
Reply by CedronMay 13, 2017

"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

[ - ]
Reply by Rick LyonsMay 13, 2017

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!