Digital Signal Processing

Lab 2: Discrete Time Systems

Convolution

Matlab provides a built-in function called conv that computes the convolution between two finite-
duration sequences x and h. The conv function assumes that the two sequences begin at n=0 and is
invoked by
y=conv(x,h);

The conv function neither provides nor accepts any timing information specifying the support of the
sequences. We define a specifically modified function called conv_m, which performs the convolution of
sequences with arbitrary support.

function [y,ny] = conv_m(x,nx,h,nh)

% Modified convolution routine for signal processing
% --------------------------------------------------
% [y,ny] = conv_m(x,nx,h,nh)
% [y,ny] = convolution result
% [x,nx] = first signal
% [h,nh] = second signal
%
nyb = nx(1)+nh(1); nye = nx(length(x)) + nh(length(h));
ny = [nyb:nye];
y = conv(x,h);

The following code is an example of explicit implementation of a function that computes the convolution
between two causal sequences x[n] and h[n]:

function y = myconv(x,h)
% y=myconv(x,h) direct computation of the convolution between x and h

L=length(x)+length(h)-1; % output length

y=zeros(1,L); % output sequence

x(L)=0; %zero padding of x till length L

h(L)=0; %zero padding of h till length L

1
 for i=1:L
for j=1:i
y(i)=y(i)+x(j)*h(i-j+1); % output computation
% Beware of the indexes (starting from 1 and not from 0 !)
end;
end;

Crosscorrelation

Given two real-valued sequences x[n] and y[n] of finite energy, the crosscorrelation of x[n] and y[n] is a
sequence rxy [n] defined as:
rxy [n] = ∑∞
k=−∞ x[k] y[k − n]
which is used to determine the similarity of x[n] and y[n] depending on a lag n. By comparing this definition
with the formula of the convolution, we find that the crosscorrelation can be also expressed as a convolution:
rxy [n] = x[n] ∗ y[−n]
Therefore the crosscorrelation can be computed by means of a convolution, if the sequences are of finite
duration.

Difference equations

An LTI discrete system can be described by a linear constant coefficient difference equation of the form
𝑁 𝑀

∑ 𝑎𝑘 𝑦[𝑛 − 𝑘] = ∑ 𝑏𝑚 𝑥[𝑛 − 𝑚] , ∀𝑛
𝑘=0 𝑚=0

A Matlab function called filter is available to solve difference equations numerically, given the input
x[n] and the difference equation coefficients. In its simplest form this function is used as follows
b = [b0, b1, ..., bM]; a = [a0, a1, ..., aN];
y = filter(b,a,x);
The output sequence y is computed on the same time interval as the input sequence.

Example: Moving-Average Filter

A moving-average filter is a common method used for smoothing noisy data. This example uses the filter
function to compute averages along a vector of data.
A moving-average filter slides a window of length M along the data, computing averages of the data
contained in each window. The following difference equation defines a moving-average filter of a vector :

2
 𝑀−1
1
𝑦[𝑛] = ∑ 𝑥[𝑛 − 𝑚]
𝑀
𝑚=0
Create a 1-by-100 row vector of sinusoidal data that is corrupted by random noise.
t = linspace(-pi,pi,100);
x = sin(t) + 0.25*rand(size(t));

For a window size of 5, compute the numerator and denominator coefficients for the rational transfer
function.
M = 5;
b = (1/M)*ones(1,M);
a = 1;

Find the moving average of the data and plot it against the original data.
y = filter(b,a,x);
plot(t,x)
hold on
plot(t,y,'r')
legend('Input Data','Filtered Data')
hold off

Difference equations with initial conditions

The filter function of Matlab assumes that the system specified by the coefficients b and a derived
from the difference equations is initially in zero-state conditions, i.e. y[n]=0 and x[n]=0 for n<0.
If this is not the case, we can take into account different initial conditions by specifying a fourth parameter
of the function filter, which represents the initial state of the filter.
The initial state of the filter is obtained by means of the Matlab function filtic.
The usage is the following (see help of filter and filtic):
ic = filtic(b, a, ypast, xpast);
where ypast and xpast are vectors of past values of y and x (for n=-1, n=-2, …)
y = filter(b, a, x, ic);

If xpast is not specified, it is assumed that x[n]=0 for n<0.

Access to final conditions fc is obtained using [y,fc]=filter(b,a,x,ic).

Example
Given the difference equation
y[n] =x[n] −0.2y[n − 1] + 0.6y[n − 2] with initial conditions y[-1]=2, y[-2]=-1
compute the output y[n] for 0≤n≤50 when the input is x[n]=0.7nu[n].

3
b= [1];
a = [1 0.2 -0.6];
n = 0:50;
x = (0.7).^n;
ic_y=[2 -1]; % initial conditions on y at time instants [-1 -2]
ic_x=[0 0]; % initial conditions on x at time instants [-1 -2]
ic=filtic(b, a, ic_y, ic_x); % set the proper filter status for n=0
y = filter(b, a, x, ic); % apply filter with initial conditions
stem(n,y)

Example Determine the zero input response of the system defined by the difference equation
𝑦[𝑛] + 1.12𝑦[𝑛 − 1] = 0.1𝑥[𝑛] + 0.2𝑥[𝑛 − 1] with initial condition y[-1] = 1.

b = [0.1, 0.2];
a = [1, 1.12];
ypast = [1]; % Initial conditions for output
ic = filtic(b,a,ypast); % Finding initial conditions for the system
yzi = filter(b,a,zeros(1,20),ic); % Zero Input response
stem(yzi)

The Matlab command impz computes the impulse response of the system corresponding to a difference
equation with coefficients a and b:
h = impz(b,a,N);
It computes the first N samples of the impulse response. If no output arguments are given, the impz
function plots the response in the current figure window using the stem function.

4
Examples
1. The convolution between the sequences x[n] and h[n] defined as follows
x = [3, 11, 7, 0, -1, 4, 2]; nx = [-3:3];
h = [2, 3, 0, -5, 2, 1]; nh = [-1:4];
is computed by
[y,ny] = conv_m(x,nx,h,nh);

2. Given the following difference equation

y[n] − y[n − 1] + 0.9y[n − 2] = x[n] +x[n − 1]; ∀n

A. Calculate and plot the impulse response h[n] at n = −20, . . . , 100.

B. Calculate and plot the unit step response s[n] at n = −20, . . . , 100.
C. Is the system specified by h[n] stable?

%A.
n = [-20:100];
a = [1, -1, 0.9];
b = [1 1];
h = impz(b,a,n);
subplot(2,1,1); stem(n,h);
title('Impulse Response'); xlabel('n'); ylabel('h(n)')
%B.
x = stepseq(0,-20,100);
s = filter(b,a,x);
subplot(2,1,2); stem(n,s)
title('Step Response'); xlabel('n'); ylabel('s(n)')

%C
%compare partial sums till n=100, 1000, 10000 …
h=impz(b,a,[-20:1000]);
sum(abs(h))
h=impz(b,a,[-20:10000]);
sum(abs(h))

3. Consider the system given by

𝑦[𝑛] − 0.5𝑦[𝑛 − 1] + 0.75𝑦[𝑛 − 2] = 2.5𝑥[𝑛] + 3𝑥[𝑛 − 1] + 2𝑥[𝑛 − 2]
and the input sequences x1[n]=cos(2πn/10) and x2[n]=sin(4πn/5). Verify the linearity of the system with
different initial conditions for y[-1] and y[-2] (e.g. {y[-1]=0 y[-2]=0} or {y[-1]=-1 y[-2]=2})

5
 n = 0:40;
% Input sequences:
x1 = cos(2*pi*0.1*n);
x2 = sin(2*pi*0.4*n);
%Weights:
C = -2;
D = 5;

x = C*x1 + D*x2;

b = [2.5 3 2];
a = [1 -0.5 0.75];
y_ic = [0 0]; % Set initial conditions for y
ic = filtic(b,a,y_ic); % Find initial conditions of the filter
y1 = filter(b,a,x1,ic); % Compute the output y1[n]
y2 = filter(b,a,x2,ic); % Compute the output y2[n]
y = filter(b,a,x,ic); % Compute the output y[n]
yt = C*y1 + D*y2;
d = y - yt; % Compute the difference of the outputs

% Plot the outputs and the difference signal

subplot(3,1,1)
stem(n,y);
ylabel('Amplitude');
title('Output when input is: C \cdot x_{1}[n] + D \cdot x_{2}[n]');
subplot(3,1,2)
stem(n,yt);
ylabel('Amplitude');
title('Weighted Output: C \cdot y_{1}[n] + D \cdot y_{2}[n]');
subplot(3,1,3)
stem(n,d);
xlabel('Time index n');
ylabel('Amplitude');
title('Difference');

4. Echo generation and reverberation: When music is performed in a concert hall, a torrent of echoes
from the various surfaces in the room strikes the ear, producing the impression of space to the listener. More
specifically, the sound reaching the listener consists of several components: direct sound, early reflections,
and reverberations. The early reflections correspond to the first few reflections off the wall, whereas the
reverberation is composed of densely packed late reflections. Music recorded in an almost anechoic studio,
using microphones placed close to the instruments, and played at home or in a car does not sound natural.
The typical solution to this problem is to create and add some amount of artificial reverberation to the
original recording before distribution.

6
A single echo is easily generated using the FIR filter:
y[n] = x[n] + ax[n − D], −1 < a < 1
where x[n] is the original signal, D is the round-trip delay in number of sampling intervals, and a is the
attenuation factor due to propagation and reflection. If the delay τ = D/Fs is greater than approximately 40
ms, an echo will be heard. A second echo will be given by a2x[n−2D], a third by a3x[n−3D], and so on.
Therefore, a multiple echo generating FIR filter is given by
y[n] = x[n] + ax[n − D] + a2x[n − 2D] + a3x[n − 3D]+· · ·
which has impulse response
h[n] = δ[n] + aδ[n − D] + a2δ[n − 2D] + a3δ[n − 3D]+· · ·
This filter generates an infinite sequence of echoes having exponentially decaying amplitudes and spaced
D sampling periods apart. A more efficient recursive implementation is
given by:
y[n] = ay[n − D] + x[n] , − 1 < a < 1
The condition −1 < a < 1 assures the stability of the system. Such simple filters provide the basic building
blocks of more sophisticated digital reverberators.

Implementation of reverberation: examples

4.1: A recursive implementation of reverberation is given by:

y[n] = x[n] + ay[n − D],
where D = τ·Fs is the delay in sampling interval given the delay τ in seconds and sampling rate Fs and a is
an attenuation factor. To generate digital reverberation we will use the sound file handel which is recorded
at Fs = 8192 samples per second.
(a) For τ = 50 ms and a = 0.7, obtain a difference equation for the digital reverberation and process the
sound in handel.
(b) Repeat (a) for τ = 100 ms.
(c) Repeat (a) for τ = 500 ms.
(d) Which implementation sounds more natural?

close all; clc

load('handel.mat')
a = 0.7; % specify attenuation factor
tau = 50e-3; % Part (a)
% tau = 100e-3; % Part (b)
%tau = 500e-3; % Part (c)
D = floor(tau*Fs); % compute delay
yd = filter(1,[1 zeros(1,D-1),-a],y);
sound(yd,Fs); % play the signal yd

7
4.2. Convolution using a recorded impulse response to simulate reverberation
% Room reverb simulation
[h,Fs] = audioread('handclap.wav'); % Recording of clap in hallway
soundsc(h,Fs)
[d,Fs] = audioread('speech.wav'); % Speech example
soundsc(d,Fs)
% Convolve them
y=(conv(d,h));
% Listen to simulated reverberation
soundsc(y,Fs)

5. The cross-correlation between two sequences can be used as a Time Delay Estimator. Let
x[n] = [3, 11, 7,0,−1, 4, 2] be a prototype sequence, and let y[n] be its noise-corrupted-and-shifted version:
y[n] = x[n – 2] + w[n]
where w[n] is a Gaussian sequence with mean 0 and variance 1. Compute the crosscorrelation between y[n]
and x[n].

x = [3, 11, 7, 0, -1, 4, 2]; nx=[-3:3]; % given signal x[n]

subplot(311); stem(nx,x);title('x[n]');
[y,ny] = sigshift(x,nx,2); % obtain x[n-2]
w = randn(1,length(y)); nw = ny; % generate w[n]
[y,ny] = sigadd(y,ny,w,nw); % obtain y[n] = x[n-2] + w[n]
subplot(312); stem(ny,y);title('y[n]');

[x1,nx1] = sigfold(x,nx); % obtain x[-n]

[rxy,nrxy] = conv_m(y,ny,x1,nx1); % crosscorrelation
subplot(3,1,3);stem(nrxy,rxy)
xlabel('lag variable l')
ylabel('rxy');title('Crosscorrelation between x and y')

We observe that the crosscorrelation indeed peaks at l= 2, which implies that y[n] is similar to x[n]
shifted by 2 samples. This approach can be used in applications like radar and sonar signal processing in
identifying and localizing targets.
Note that the signal-processing toolbox in MATLAB also provides a function called xcorr for sequence
correlation computations. In its simplest form:
>> xcorr(x,y)
it computes the crosscorrelation between vectors x and y, while
>> xcorr(x)
computes the autocorrelation of vector x. It generates results that are identical to the one obtained from the
proper use of the conv_m function. However, the xcorr function cannot provide the timing (or lag)
information (as done by the conv_m function), which then must be obtained by some other means.

8
6. Compute and plot the output of the discrete-time system. y[n] = 5y[n − 1] + x[n], y[−1] = 0, for x[n] =
u[n], 0 ≤ n ≤ 100. Based on these results can you make a statement regarding the stability of the system?
Hint: What happens when n increases?

close all; clc

n = 0:100;
x = ones(1,length(n));
y = filter(1,[1 -5],x);
% Plot:
stem(n,y,'fill')
xlabel('n');
title('y[n] = 5y[n-1]+x[n], y[-1]=0')

7. Consider the system y[n] = y[n − 1] + y[n − 2] + x[n], y[−1] = y[−2] = 0.

(a) Compute and plot the impulse response, for 0 ≤ n ≤ 100, using the function filter.
(b) Can you draw any conclusions about the stability of this system from the results in (a)?

close all; clc

% Part (a):
n = 0:100;
[x n]=impseq(n(1),n(1),n(end));
y = filter(1,[1 -1 -1],x);
% Plot:
stem(n,y, 'fill')
xlabel('n'); ylabel('y[n] ')
title('LCCDE: y[n] = y[n-1] + y[n-2] + x[n], y[-1] = y[-2] = 0')

8. Given the system specified by

y[n] = 1.15y[n − 1] − 1.5y[n − 2] + 0.7y[n − 3] − 0.25y[n − 4] + 0.18x[n] + 0.1x[n − 1] + 0.3x[n − 2] +
0.1x[n − 3] + 0.18x[n − 4].

(a) Compute and plot the impulse response h[n], 0 ≤ n ≤ 100 using the function h=impz(b,a,N).
(b) Compute and plot the output y[n], if x[n] = u[n], 0 ≤ n ≤ 100 using the function y=filter(b,a,x).
(c) Compute and plot the output y[n], if x[n] = u[n], 0 ≤ n ≤ 100 using the function y=conv(h,x).
(d) Compute and plot the output y[n], if x[n] = u[n], 0 ≤ n ≤ 100 using the function y=filter(h,1,x).

close all; clc

n = 0:100;
b = [0.18 0.1 0.3 0.1 0.18];
a = [1 -1.15 1.5 -0.7 0.25];
% Part (a):
h = impz(b,a,length(n));
% Part (b):
[u n]=stepseq(n(1),n(1),n(end));
y1 = filter(b,a,u);

9
% Part (c):
y2 = conv(h,u);
% Part (d):
y3 = filter(h,1,u);
%Plot
subplot(221)
stem(n,h, 'fill')
xlabel('n'); ylabel('h[n] ');
title('Impulse Response h[n] ');
subplot(222)
stem(n,y1, 'fill')
xlabel('n'); ylabel('y_1[n] ');
title('Unit Step Response: filter(b,a,x)');
subplot(223)
stem(n,y2(1:length(n)),'fill')
xlabel('n'); ylabel('y_2[n]');
title('Unit Step Response: conv(h,x)');
subplot(224)
stem(n,y3,'fill')
xlabel('n'); ylabel('y_3[n]');
title('Unit Step Response: filter(h,1,x)');

9. Consider the finite duration sequences x[n] = u[n] − u[n − N] and h[n] = n(u[n] − u[n − M]), M ≤ N.
Assume N = 10 and M = 5 and compute the convolution y[n] between x[n] and h[n].

N = 10; M = 5;
n = 0:N-1;
x=ones(1,length(n));
h=n(1:M).*ones(1,M);
[y ny] = conv_m(h,0:M-1,x,n);
% Plot:
stem(ny,y, 'fill')
axis([ny(1)-1,ny(end)+1,min(y)-1,max(y)+1])
xlabel('n'); ylabel('y[n]');
title('y[n] = h[n]*x[n]')

10. A system is described by the difference equation y[n] = x[n] − 0.9y[n − 1], x[n] = u[n] – u[n – 10]
What is the output of the system? (hint: plot y[n] for a suitable range of indexes n)
b = [1]; a = [1,0.9];
n = -5:50; x = stepseq(0,-5,50) - stepseq(10,-5,50);
y = filter(b,a,x);
stem(n,y); title('Output sequence')
xlabel('n'); ylabel('y(n)'); axis([-5,50,-0.5,1])

10