Forums

Simulation of DTFT in MATLAB

Started by krishna_sun82 April 6, 2006
How do I have to simulate the Discrete time Fourier Transform in MATLAB? Do
I have to make use of freqz function for it? Thanks in advance.

- Krishna
krishna_sun82 wrote:
> How do I have to simulate the Discrete time Fourier Transform in MATLAB? Do > I have to make use of freqz function for it? Thanks in advance.
For finite sequences, the DTFT is just the z-transform evaluated on the unit circle. So yes, if x is your vector then you get the DTFT evaluated at the frequencies F (in Hz) by H = freqz(x,1,F,Fs), where Fs is the sampling frequency.
On Thu, 06 Apr 2006 00:53:27 -0500, "krishna_sun82"
<krishna.sun@gmail.com> wrote:

>How do I have to simulate the Discrete time Fourier Transform in MATLAB? Do >I have to make use of freqz function for it? Thanks in advance. > >- Krishna
You simulate the DTFT in MATLAB by zero padding a time-domain sequence with many many zeros. Then perform the DFT on that extended-length sequence by whatever DFT method [fft() or freqz()] you wish. [-Rick-]
Rick Lyons wrote:
> On Thu, 06 Apr 2006 00:53:27 -0500, "krishna_sun82" > <krishna.sun@gmail.com> wrote: > > >How do I have to simulate the Discrete time Fourier Transform in MATLAB? Do > >I have to make use of freqz function for it? Thanks in advance. > > > >- Krishna > > You simulate the DTFT in MATLAB by zero padding > a time-domain sequence with many many zeros. > > Then perform the DFT on that extended-length sequence > by whatever DFT method [fft() or freqz()] you wish.
he might not need to zero-pad if he can define the whole extended length sequence meaningfully. but the point is the same. since the output of the DTFT is a continuous-frequency function with infinite number of points (usually between w=-pi and w=+pi), and since computers have trouble dealing with any infinite quantity of any object, if you're gonna simulate, you have to decide on some finite number that is close enough. that would be the DFT "N". the other thing i would recommend is to fftshift() the input to the fft() (to avoid that phase being off by +/- pi for odd X[k]) and to also fftshift() the output of the fft() so that the negative frequencies components come first. (if MATLAB had the ability to do negative indices, a version of the DFT could be written so that none of this swapping of the two halves would be necessary.) N = some_really_big_number_likely_a_power_of_two; omega = linspace(-pi, pi(1-2/N), N); x = some_really_big_array_of_size_N; X = fftshift( fft( fftshift(x) ) ); plot(omega, X); % this is N points of the DFT in a linear scale r b-j r b-j
robert bristow-johnson wrote:
...
> but the point is the same. since the output of the DTFT is a > continuous-frequency function with infinite number of points (usually > between w=-pi and w=+pi), and since computers have trouble dealing with > any infinite quantity of any object, if you're gonna simulate, you have > to decide on some finite number that is close enough. that would be > the DFT "N".
freqz() can evalute the DTFT on any arbitrary frequency points (not necessarily evenly spaced, as is the result of the FFT).
Andor wrote:
> > freqz() can evalute the DTFT on any arbitrary frequency points (not > necessarily evenly spaced, as is the result of the FFT).
that's true, but it would still have to be evalutated at a finite selection of points, though they need not be uniformly spaced. i made a small screwup in the MATLAB pseudo-code in my previous post: N = some_really_big_number_likely_a_power_of_two; omega = linspace(-pi, pi(1-2/N), N); x = some_really_big_array_of_size_N; X = fftshift( fft( fftshift(x) ) ); plot(omega, abs(X)); % this is N points of the DFT magnitude in a linear scale pause: plot(omega, angle(X)); % this is N points of the DFT wrapped phase in a linear scale i *still* recommend the two fftshift() operations. there is a way to do split screen plots, but i forgot it. r b-j
Thank you all for your replies. I tried using both freqz and fftshift. In
my case, there is not much difference in the end result. I don't really
use arbitrary frequency values, but still freqz provides that flexibility
over fftshift.

>i *still* recommend the two fftshift() operations. >there is a way to do split screen plots, but i forgot it. > >r b-j >
Robert, can you please tell me why and in which case you would prefer fftshift over freqz? Is there any mathematical reason that I fail to comprehend? For split screen plots, you may use subplot(m, n, p) function. One more doubt: is there any way to find the inverse of the DTFT in MATLAB? Also, DTFT output is actually an infinite sequence, but in a computer we can get only finite points. To calculate the inverse DTFT, we have to integrate the scaled DTFT output over 2pi. But the computer output for DTFT contains finite points. If there is no instant function in MATLAB, can I have summation in the place of integration? What other factors do I need to consider to have a reasonably accurate output?
having trouble sleeping, so...

krishna_sun82 wrote:
> Thank you all for your replies. I tried using both freqz and fftshift. In > my case, there is not much difference in the end result. I don't really > use arbitrary frequency values, but still freqz provides that flexibility > over fftshift. > > Robert, can you please tell me why and in which case you would prefer > fftshift over freqz? Is there any mathematical reason that I fail to > comprehend?
it's not just fftshift() vs. freqz(). it's fftshift( fft( fftshift() ) ) vs. freqz(). all fftshift() does is swap the first half with the second half of the array. suppose your chunk of data of which spectrum you are interested in was yanked out of a larger stream of samples (that is essentially applying a window function, even if it's the rectangular window function). also, suppose you used a "better" window function, like a Hamming or Kaiser, so that the samples in the middle have more amplitude than the samples at the two ends. now, where would you like your "y-axis" to go through (where t=0 is), the middle of that chunk or one of the ends? that's the reason to use fftshift on the data going *in* to the FFT, so that the data in the second half of the input is treated as the data left of the t=0 mark, the negative time data. doing it on the output is so that the negative frequency output comes before, to the left of, the positive frequency output.
> > One more doubt: is there any way to find the inverse of the DTFT in > MATLAB?
the inverse DTFT is an integral (of a continuous-frequency function). unless you're doing it mathematically or symbolically, evaluating an integral of a continuous function would entail an uncountably infinite number of additions, something pretty hard to make a computer do. so like the DFT is an evaluation of the DTFT at a finite number of equally spaced points, the iDFT is an approximation of the iDTFT where the integral is replaced by a Riemann summation. so if you have a continuous-frequency function of w, X(e^jw) that is the DTFT of an unknown discrete-time sequence, x[n], and you want to numerically compute some approximation to x[n], then uniformly sample X(e^jw) a gazillion (N) times from w=-pi to w=+pi, swap the halfs (using fftshift) and inverse DFT that result and you'll have a gazillion samples of x[n]. but since the real x[n] might never have been periodic (which is why it is describable with the DTFT of it), the output of the iDFT will be one period of a periodic sequence (of period N) that should very well approximate your aperiodic x[n] for small enough n. that is the consequence of sampling X(e^jw), since X(e^jw) has more information in it than the samples of it: X[k] = X(e^jw) where w = 2*pi*k/N .
> Also, DTFT output is actually an infinite sequence,
no, that's the INPUT to the DTFT. the output of it is a continuous-frequency function from -pi to +pi.
> but in a computer we > can get only finite points. To calculate the inverse DTFT, we have to > integrate the scaled DTFT output over 2pi. But the computer output for > DTFT contains finite points.
geez, i wish i would read through your post before typing. it might save me some keystrokes.
> If there is no instant function in MATLAB, > can I have summation in the place of integration?
yes, you have to. but what you're summing are samples of the DTFT output (with something multipying them) which has the effect of periodically repeating and overlap-adding the original infinite sequence x[n] resulting in a periodic (or aliased) extension of x[n].
> What other factors do I > need to consider to have a reasonably accurate output?
well your original x[n] that is an infinite sequence must decay to zero rapidly enough in both directions (call those the "tails" of x[n]) so that the DTFT summation converges to finite values of X(e^jw). when you grab a finite segment of x[n] to pass to a DFT (which is what you are using to approximate the DTFT), you are windowing it. to minimize the effect of windowing, make sure the window is wide enough to include most of the tails (so that what gets cut off is very small in amplitude). perhaps, just read the enormous back-and-forth we've been having in the "Periodicity of the DFT" subthread we've had very recently. i don't want to repeat myself all over again. r b-j