518 Practical MATLAB® Applications for Engineers
Example 5.6
Let the input sequence to a system be given by f(n) = [1 0 3 2], and its impulse system
response be h(n) = [−1 0 3 4 6 3].
Create the script fi le convolutions that returns
a. The system output g(n)= f(n) ⊗ h(n), using the direct time convolution, as well as
the DFT of g(n)
b. The system output g(n) by evaluating the IDFT of the product of the DFT of f(n) with
the DFT of h(n)
c. The estimate of the error of the methods used in part a with respect to part b by
means of a plot
ANALYTICAL Solution
Note that DFT involves circular convolution, therefore the sequences f(n) and h(n) must
be augmented by appending zeros to conform to a length of nine (6 + 4 − 1 = 9) for
each sequence.
MATLAB Solution
% Script file: convolutions
fn = [1 0 3 2];
hn = [-1 0 3 4 6 3];
gn = conv(fn,hn); % direct convolution
n = 0:8;
fn _ aug = [fn zeros(1,5)];
hn _ aug = [hn zeros(1,3)];
DTF _ fn _ aug = fft(fn _ aug); % DFTs by zero padding
DTF_hn_aug = fft(hn_aug);
DTF _ aug = DTF _ fn _ aug.*DTF _ hn _ aug;
ggn = ifft(DTF _ aug); % linear convolution
figure(1)
subplot (3,1,1)
stem (n,gn); hold on; plot(n,gn); ylabel (‘Amplitude’);
xlabel(‘time index n’)
title(‘ Linear convolution [g(n)=conv(f(n),h(n))] vs. n’)
subplot(3,1,2)
DFT _ lin = fft(gn,9);
stem (n,abs(DFT _ lin));hold on; plot(n,abs(DFT _ lin));
xlabel (‘frequency W’);
ylabel(‘Magnitude’);
title (‘Magnitude of DFT[g(n)] vs.W’)
subplot(3,1,3)
stem (n, angle(DFT _ lin)); hold on; plot(n, angle(DFT _ lin));
title (‘Phase of DFT[g(n)] vs. W’); xlabel(‘frequency W’)
ylabel(‘Angle in rad.’)
figure(2)
subplot(3,1,1)