DTFT of a signal

Started by 6 years ago9 replieslatest reply 6 years ago1150 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); gridxlabel ('frequency in pi units'); ylabel('|X|')title ('Magnitude part');subplot (2,1,2);plot (w/pi, angX/pi); gridxlabel ('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:

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

[ - ]