PRACTICAL MATLAB® FOR ENGINEERS PRACTICAL MATLAB

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)