Computing the DTFT of a signal in Matlab depends on
a) if the signal is finite duration or infinite duration
b) do we want the numerical computation of the DTFT or a closed form expression.
In the examples that follow, u[n] is the discrete time unit step function, i.e.,
u[n] = 1, n >= 0
u[n] = 0, n < 0
Case 1: finite duration, numerical computation
Let:
x1[n] = (0.5^n)*(u[n] - u[n - 8])
x2[n] = (1.5^n)*(u[-n] - u[-n + 5])
x[n] = x1[n] + x2[n];
x[n] is finite duration, i.e., x[n] = 0 for n < -4 and n > 7.
u = @(n) double(n>=0); % discrete unit step
x1 = @(n) (0.5.^n).*(u(n) - u(n-8));
x2 = @(n) (1.5.^n).*(u(-n) - u(-n - 5));
x = @(n) x1(n) + x2(n);
Numerical computation of the DTFT can be accomplished in a couple of ways. Perhaps the easiest is to use the Signal Processing Toolbox function freqz. However, using freqz in this way assumes that the first element of the first input to freqz corresponds to n = 0. But, in this example the first nonzero element of x[n] is at n = -4. We have to use the shift property of the DTFT after calling freqz to ensure we get the correct phase.
Let y[n] be a version of x[n] shifted to the right by m samples such that all of the non-zero elements of y[n] are at times n >= 0.
y[n] = x[n - m], m > 0
In our case, we have to shift x[n] to the right by at least 4 samples.
The DTFT of y[n] is related to the DTFT of x[n] as
Y(w) = X(w)*exp(-1j*w*m) (w in rad)
Therefore, we can use freqz to compute Y(w), then invert the above to compute X(w)
X(w) = Y(w)*exp(1j*w*m)
Compute x[n] over its duration
n = -4:7;
xval = x(n);
Compute freqz of xval, which are really the values of yval as far as freqz is concerned. Use the 'whole' option to evaluate around the entire unit circle.
[h,w] = freqz(xval,1,2048,'whole');
Phase adjustment for the time domain shift implicit to freqz
h = h.*exp(1j*w*4);
h is basically one period of the DTFT of x[n]. Plot it
figure
plot(subplot(211),w,abs(h)); grid, ylabel('Magnitude'),ylim([0 5]),xlim([0 2*pi])
plot(subplot(212),w,angle(h)); grid, ylabel('Phase (rad)'), xlabel('w (rad)'),xlim([0 2*pi])
Another option, which might be slightly more numerically robust, is to use fft with lots of zero padding to compute one period of the DTFT at lots of frequencies. fft assumes the first point in the sequence corresponds to n = 0, so we have to circshift after zero padding and then apply fft
h = fft(circshift([xval zeros(1,2048-numel(xval))],-4),2048);
w = (0:2047)/2048*2*pi;
figure
plot(subplot(211),w,abs(h)); grid, ylabel('Magnitude'),ylim([0 5]),xlim([0 2*pi])
plot(subplot(212),w,angle(h)); grid, ylabel('Phase (rad)'), xlabel('w (rad)'), xlim([0 2*pi])
Another option would be write a function to compute the DTFT sum directly from its definition at frequencies of interest.
Case 2: Finite duration, closed form expression.
This case is straightforward using the Symbolic Math Toolbox
clear all
syms n integer
u(n) = kroneckerDelta(n)/2 + heaviside(n); % discrete time unit step using default sympref
x1(n) = (sym(0.5)^n)*(u(n) - u(n-8));
x2(n) = (sym(1.5)^n)*(u(-n) - u(-n - 5));
x(n) = x1(n) + x2(n);
We know that x[n] is zero outside the interval -4 <= n <= 7, so we can use the DTFT sum over that interval explicitly
syms w
nval = -4:7;
X(w) = sum(x(nval).*exp(-1j*w*nval))
Plot one period
figure
fplot(subplot(211),abs(X(w)),[0 2*pi]), grid, ylabel('Magnitude'), ylim([0 5])
fplot(subplot(212),angle(X(w)),[0 2*pi]), grid, ylabel('Phase (rad)'), xlabel('w (rad)')
Case 3: Infinite duration, closed form expression
Let x[n] be similar to above, but instead of cutting x[n] off at n = -5 and n = 8, we'll let it run out to abs(n) = inf.
clear all
syms n integer
u(n) = kroneckerDelta(n)/2 + heaviside(n); % discrete time unit step using default sympref
x1(n) = (sym(0.5)^n)*(u(n));
x2(n) = (sym(1.5)^n)*(u(-n));
x(n) = x1(n) + x2(n);
Matlab's Symbolic Math Toolbox doesn't offer a symbolic DTFT function, but it does offer the unilateral z-transform via ztrans/iztrans. We can use ztrans to compute the bilateral z-transform of x1[n], which is the same as its unilateral transform because x1[n] is causal. We can use ztrans to to compute the bilateral z-transform of x2[n] using the z-transform time reversal property because x2[n] is strictly non-causal. Because the bilateral z-transforms of x1[n] and x2[n] both have a region of convergence that includes the unit circle, so does their sum, in which case the DTFT of x[n] can be computed by evaluating its bilateral z-tansform around the unit circle.
syms z
X1(z) = ztrans(x1(n),n,z); % z-transform of causal signal, ROC: abs(z) > 0.5
X2(z) = simplify(subs(ztrans(x2(-n),n,z),z,1/z)); % z-transform of non-causal signal, ROC: abs(z) < 1.5
X(z) = X1(z) + X2(z); % ROC: 0.5 < abs(z) < 1.5
syms w
X(w) = X(exp(1j*w)); % DTFT of x[n]
Plot one period
figure
fplot(subplot(211),abs(X(w)),[0 2*pi]), grid, ylabel('Magnitude'), ylim([0 5])
fplot(subplot(212),angle(X(w)),[0 2*pi]), grid, ylabel('Phase (rad)'), xlabel('w (rad)'), ylim([-0.2 0.2])
Case 4: Infinite duration, numerical computation
Because x[n] is infinite duration we can't really do a numerical computation, which would need an infinite number of samples. But x[n] decays to zero in both directions, so we can pick a duration of nlower <= n <= nupper, where for values of n outside those limits it's safe to assume that x[n] can be approximated as x[n] = 0.
We'll use upper and lower bounds of 20
nval = -20:20;
xval = double(x(nval));
Now we can use any of the numerical techniques from Case 1 with xval being elements of the finite duration approximation to x[n]
[hval,wval] = freqz(xval,1,2048,'whole');
hval = hval.*exp(1j*wval*20);
Plot one period
figure
plot(subplot(211),wval,abs(hval)); grid, ylabel('Magnitude'),ylim([0 5]),xlim([0 2*pi])
plot(subplot(212),wval,angle(hval)); grid, ylabel('Phase (rad)'), xlabel('w (rad)'), xlim([0 2*pi]), ylim([-0.2 0.2])
All of this becomes much simpler if the time domain signal, x[n], is causal, i.e., 0 for n < 0.
