Discrete Signal Processing (DSP):

Linear Time-Invariant (LTI) Systems

By Prof. Seungchul Lee
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH

Table of Contents

1. Linear Time-Invariant (LTI) Systems

A discrete-time system ${H}$ is a transformation (a rule or formula) that maps a discrete-time input signal $x$ into a discrete-time output signal $y$

System examples

$$ \begin{align*} &\text{Identity} \quad &y[n] &= x[n] \quad &\forall n\\ &\text{Scaling} \quad &y[n] &= 2\,x[n] \quad &\forall n\\ &\text{Offset} \quad &y[n] &= x[n]\, + \,2 \quad &\forall n\\ &\text{Square signal} \quad &y[n] &= (x[n])^2 \quad &\forall n\\ &\text{Shift} \quad &y[n] &= x[n + 2] \quad &\forall n\\ &\text{Decimate} \quad &y[n] &= x[2n] \quad &\forall n\\ &\text{Square time} \quad &y[n] &= x[n^2] \quad &\forall n \end{align*} $$

1.1. Linear Systems

A system ${H}$ is linear if it satisfies the following two properties:

1) Scaling

$${H} \{\alpha x \} = \alpha {H}\{x\} \quad \forall \, \alpha \in \mathbb{C}$$

2) Additivity

$$ \text{If}\,\, y_1 = {H} \{x_1\} \,\, \text{and}\,\, y_2 = {H} \{ x_2 \} \, \, \text{then } {H} \{x_1 + x_2\} = y_1 + y_2$$

1.1.1 Matrix Multiplication and Linear Systems

• Matrix multiplication (aka Linear Combination) is a fundamental signal processing system

• Matrix multiplications are linear systems

$$\begin{align*} y &= Hx \\\\ y[n] &= \sum_{m} [H]_{n,m} x[m] = \sum_{m}h_{n,m} x[m] \end{align*}$$

$\quad\;$ where $h_{n,m} = [H]_{n,m}$ represents the row-$n$, column-$m$ entry of the matrix $H$

• All linear systems can be expressed as matrix multiplications

1.1.2. Matrix Multiplication and Linear Systems in Pictures

• System output as a linear combination of columns

• System output as a sequence of inner products

1.2. Time-Invariant Systems

For infinite-length signals

• A system ${H}$ processing infinite-length signals is time-invariant (shift-invariant) if a time shift of the input signal creates a corresponding time shift in the output signal

For finite-length signals

• A system ${H}$ processing infinite-length signals is time-invariant (shift-invariant) if a circular time shift of the input signal creates a corresponding circular time shift in the output signal

1.3. Linear Time-Invariant (LTI) Systems

A system $H$ is linear time-invariant (LTI) if it is both linear and time-invariant

We will only consider Linear Time-Invariant (LTI) systems.

$$\begin{align*} \text{Identity} &\qquad y[n] = x[n] \quad &\forall n &\qquad \text{Linear} &\qquad \text{Time Invariant}\\ \text{Scaling} &\qquad y[n] = 2\,x[n] \quad &\forall n &\qquad \text{Linear} &\qquad \text{Time Invariant}\\ \text{Offset} &\qquad y[n] = x[n]\, + \,2 \quad &\forall n &\qquad \text{Non Linear} &\qquad \text{Time Invariant}\\ \text{Square signal} &\qquad y[n] = (x[n])^2 \quad &\forall n &\qquad \text{Non Linear} &\qquad \text{Time Invariant}\\ \text{Shift} &\qquad y[n] = x[n + 2] \quad &\forall n &\qquad \text{Linear} &\qquad \text{Time Invariant}\\ \text{Decimate} &\qquad y[n] = x[2n] \quad &\forall n &\qquad \text{Linear} &\qquad \text{Time Variant}\\ \text{Square time} &\qquad y[n] = x[n^2] \quad &\forall n &\qquad \text{Linear} &\qquad \text{Time Variant}\\ \end{align*}$$

1.3.1. Matrix Multiplication and LTI Systems (Infinite-Length Signals)

• When the linear system is also shift invariant, $H$ has a special structure

• Linear system for infinite-length signals can be expressed as

$$y[n]=H\{x[n]\} = \sum_{m=-\infty}^{\infty} h_{n,m}x[m], \quad -\infty < n < \infty$$

• Enforcing time invariance implies that for all $q \in \mathbb{Z}$

$$H\{x[n-q]\} = \sum_{m=-\infty}^{\infty} h_{n,m}x[m-q] = y[n-q]$$

• Change of variables: $n'=n-q$ and $m'=m-q$

$$H\{x[n']\} = \sum_{m=-\infty}^{\infty} h_{n'+q,m'+q}x[m'] = y[n']$$

• Comparing first and third equations, we see that for an LTI system

$$h_{n,m} = h_{n+q,m+q}$$

$$H = \begin{bmatrix}&\vdots&\vdots&\vdots&\\ \cdots&h_{-1,-1}&h_{-1,0}&h_{-1,1}&\cdots\\ \cdots&h_{0,-1}&h_{0,0}&h_{0,1}&\cdots\\ \cdots&h_{1,-1}&h_{1,0}&h_{1,1}&\cdots\\ &\vdots&\vdots&\vdots&\\ \end{bmatrix}= \begin{bmatrix}&\vdots&\vdots&\vdots&\\ \cdots&h_{0,0}&h_{-1,0}&h_{-2,0}&\cdots\\ \cdots&h_{1,0}&h_{0,0}&h_{-1,0}&\cdots\\ \cdots&h_{2,0}&h_{1,0}&h_{0,0}&\cdots\\ &\vdots&\vdots&\vdots&\\ \end{bmatrix}$$

• Entries on the matrix diagonals are the same - Toeplitz matrix

• All of the entries in a Toeplitz matrix can be expressed in terms of the entries of the
  • 0-th column $$h[n] = h_{n,0}$$
  • Time-reversed 0-th row $$h[m] = h_{0,-m}$$

$$\begin{bmatrix}&\vdots&\vdots&\vdots&\\ \cdots&h_{0,0}&h_{-1,0}&h_{-2,0}&\cdots\\ \cdots&h_{1,0}&h_{0,0}&h_{-1,0}&\cdots\\ \cdots&h_{2,0}&h_{1,0}&h_{0,0}&\cdots\\ &\vdots&\vdots&\vdots&\\ \end{bmatrix}= \begin{bmatrix}&\vdots&\vdots&\vdots&\\ \cdots&h[0]&h[-1]&h[-2]&\cdots\\ \cdots&h[1]&h[0]&h[-1]&\cdots\\ \cdots&h[2]&h[1]&h[0]&\cdots\\ &\vdots&\vdots&\vdots&\\ \end{bmatrix} $$

• Row-$n$, column-$m$ entry of the matrix $$[H]_{n,m} = h_{n,m} = h[n-m]$$

• Note the diagonals !

1.3.2. Matrix Multiplication and LTI Systems (Finite-Length Signals)

• Linear system for signals of length $N$ can be expressed as

$$y[n]=H\{x[n]\} = \sum_{m=0}^{N} h_{n,m}x[m], \quad 0 \leq n \leq N-1$$

• Enforcing time invariance implies that for all $q \in \mathbb{Z}$

$$H\{x[(n-q)_N]\} = \sum_{m=0}^{N-1} h_{n,m}x[(m-q)_N] = y[(n-q)_N]$$

• Change of variables: $n'=n-q$ and $m'=m-q$

$$H\{x[(n')_N]\} = \sum_{m=-q}^{N-1-q} h_{(n'+q)_N,(m'+q)_N}x[(m')_N] = y[(n')_N]$$

• Comparing first and third equations, we see that for an LTI system

$$h_{n,m} = h_{(n+q)_N,(m+q)_N}$$

$$H = \begin{bmatrix}h_{0,0}&h_{0,1}&h_{0,2}&\cdots&h_{0,N-1}\\ h_{1,0}&h_{1,1}&h_{1,2}&\cdots&h_{1,N-1}\\ h_{2,0}&h_{2,1}&h_{2,2}&\cdots&h_{2,N-1}\\ \vdots&\vdots&\vdots&&\vdots\\ h_{N-1,0}&h_{N-1,1}&h_{N-1,2}&\cdots&h_{N-1,N-1}\\ \end{bmatrix}= \begin{bmatrix}h_{0,0}&h_{N-1,0}&h_{N-2,0}&\cdots&h_{1,0}\\ h_{1,0}&h_{0,0}&h_{N-1,0}&\cdots&h_{2,0}\\ h_{2,0}&h_{1,0}&h_{0,0}&\cdots&h_{3,0}\\ \vdots&\vdots&\vdots&&\vdots\\ h_{N-1,0}&h_{N-2,0}&h_{N-3,0}&\cdots&h_{0,0}\\ \end{bmatrix}$$

• Entries on the matrix diagonals are the same + circular wraparound - Circulent matrix

• All of the entries in a circulent matrix can be expressed in terms of the entries of the
  • 0-th column $$h[n] = h_{n,0}$$
  • Circularly time-reversed 0-th row $$h[m] = h_{0,(-m)_N}$$

$$\begin{bmatrix}h_{0,0}&h_{N-1,0}&h_{N-2,0}&\cdots&h_{1,0}\\ h_{1,0}&h_{0,0}&h_{N-1,0}&\cdots&h_{2,0}\\ h_{2,0}&h_{1,0}&h_{0,0}&\cdots&h_{3,0}\\ \vdots&\vdots&\vdots&&\vdots\\ h_{N-1,0}&h_{N-2,0}&h_{N-3,0}&\cdots&h_{0,0}\\ \end{bmatrix}= \begin{bmatrix}h[0]&h[N-1]&h[N-2]&\cdots&h[1]\\ h[1]&h[0]&h[N-1]&\cdots&h[2]\\ h[2]&h[1]&h[0]&\cdots&h[3]\\ \vdots&\vdots&\vdots&&\vdots\\ h[N-1]&h[N-2]&h[N-3]&\cdots&h[0]\\ \end{bmatrix}$$

• Row-$n$, column-$m$ entry of the matrix $$[H]_{n,m} = h_{n,m} = h[(n-m)_N]$$

• Note the diagonals and circulent shifts !

2. Impulse Response

• The 0-th column of the matrix $H$ has a special interpretation

• Compute the output when the input is a delta function (impulse)

$$\text{(impulse)}\;\; \delta[n] = \begin{cases} 1 \quad n = 0\\ 0 \quad \text{otherwise}\\ \end{cases}$$

• Illustrate it with an infinite-length signal

• This suggests that we call $h$ the impulse response of the system

• Output of system to delta function (impulse) is $h$. So, We call $h$ the impulse response of the system

• From $h$, we can build matrix $H$
  • Columns/rows of $H$ are the shifted versions of the 0-th column/row
  • $h$ contains all the information of the LTI system

• LTI systems are Toeplitz matrices (infinite-length signals)
  • Entries on the matrix diagonals are the same

$$y[n] = \sum_{-\infty}^{\infty}h[n-m]x[m]$$

• Let the input $x[n] = \delta[n]$ and compute

$$\sum_{-\infty}^{\infty}h[n-m]\delta[m] = h[n]$$

• The impulse response characterizes an LTI system (that is, carries all of the information contained in matrix $H$)

2.1. Examples (Infinite-Length Signals)

$$H= \begin{bmatrix}&\vdots&\color{red}{\vdots}&\vdots&\\ \cdots&h[0]&\color{red}{h[-1]}&h[-2]&\cdots\\ \cdots&h[1]&\color{red}{h[0]}&h[-1]&\cdots\\ \cdots&h[2]&\color{red}{h[1]}&h[0]&\cdots\\ &\vdots&\color{red}{\vdots}&\vdots&\\ \end{bmatrix} $$

• Impulse response of the scaling system

$$ \begin{align*} h[n] &= 2\delta[n]\\\\ H_{n,m} &= h[n-m] \end{align*} $$

%plot -s 560,200

N = 8;
nh = [0:2*N-1];
nx = [-N:N-1];

x = zeros(size(nx))';
x(N+1) = 1;

h = zeros(size(nh))';
h(1) = 2;

subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,2])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,2])

%plot -s 560,420

hc = h;
hr = [h(1),zeros(1,length(x)-1)];
H = toeplitz(hc,hr);

imagesc(H); colormap('jet'); axis('square'); axis off

• Impulse response of the shift system

$$ \begin{align*} h[n] &= \delta[n-3]\\\\ H_{n,m} &= h[n-m] \end{align*} $$

%plot -s 560,200

N = 8;
nh = [0:2*N-1];
nx = [-N:N-1];

x = zeros(size(nx))';
x(N+1) = 1;

h = zeros(size(nh))';
h(4) = 1;

subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])

%plot -s 560,420

hc = h;
hr = [h(1),zeros(1,length(x)-1)];
H = toeplitz(hc,hr);

imagesc(H); colormap('jet'); axis('square'); axis off

• Impulse response of the moving average system

$$ \begin{align*} h[n] &= \frac{1}{3}(\delta[n] + \delta[n-1] + \delta[n-2])\\\\ H_{n,m} = h[n-m] &= \frac{1}{3}(\delta[n-m] + \delta[n-m-1] + \delta[n-m-2]) \end{align*} $$

%plot -s 560,200

N = 8;
nh = [0:2*N-1];
nx = [-N:N-1];

x = zeros(size(nx))';
x(N+1) = 1;

h = zeros(size(nh))';
h(1) = 1/3;
h(2) = 1/3;
h(3) = 1/3;

subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])

%plot -s 560,420

hc = h;
hr = [h(1),zeros(1,length(x)-1)];
H = toeplitz(hc,hr);

imagesc(H); colormap('jet'); axis('square'); axis off

• Impulse response of the recursive average system

$$y[n] = x[n] + \alpha y[n-1]$$

$$ \begin{align*} h[n] &= \alpha^n u[n]\\\\ H_{n,m} = h[n-m] &= \alpha^{n-m}u[n-m] \end{align*} $$

%plot -s 560,200

N = 8;
nh = [0:2*N-1];
nx = [-N:N-1];

x = zeros(size(nx))';
x(N+1) = 1;

h = 0.8.^nh;
h = h';

subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])

%plot -s 560,420

hc = h;
hr = [h(1),zeros(1,length(x)-1)];
H = toeplitz(hc,hr);

imagesc(H); colormap('jet'); axis('square'); axis off

2.2. Examples (Finite-Length Signals)

• Entries on the matrix diagonals are the same + circular wraparound

$$h_{n,m} = h_{(n+q)_N (m+q)_N} \qquad \forall q \in \mathbb{Z}$$

$$ H = \begin{bmatrix}\color{red}{h[0]}&h[N-1]&h[N-2]&\cdots&h[1]\\ \color{red}{h[1]}&h[0]&h[N-1]&\cdots&h[2]\\ \color{red}{h[2]}&h[1]&h[0]&\cdots&h[3]\\ \vdots&\vdots&\vdots&&\vdots\\ \color{red}{h[N-1]}&h[N-2]&h[N-3]&\cdots&h[0]\\ \end{bmatrix} $$

• Impulse response of the shift system

$$ \begin{align*} h[n] &= \delta[n-2]\\\\ H_{n,m} &= h[(n-m)_N] \end{align*} $$

%plot -s 560,200

N = 8;
nh = [0:N-1];
nx = [0:N-1];

x = zeros(size(nx))';
x(1) = 1;

h = zeros(size(nh))';
h(3) = 1;

subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])

%plot -s 560,420

H = zeros(N,N);

for k = 0:N-1
    H(:,k+1) = circshift(h,k);
end

imagesc(H); colormap('jet'); axis('square'); axis off

• Impulse response of the moving average system

$$ \begin{align*} h[n] &= \frac{1}{3}(\delta[n] + \delta[n-1] + \delta[n-2])\\\\ H_{n,m} &= h[(n-m)_N] \end{align*} $$

%plot -s 560,200

N = 8;
nh = [0:N-1];
nx = [0:N-1];

x = zeros(size(nx))';
x(1) = 1;

h = zeros(size(nh))';
h(1) = 1/3;
h(2) = 1/3;
h(3) = 1/3;

subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])

%plot -s 560,420

H = zeros(N,N);

for k = 0:N-1
    H(:,k+1) = circshift(h,k);
end

imagesc(H); colormap('jet'); axis('square'); axis off

• Impulse response of the recursive average system

$$y[n] = x[n] + \alpha y[n-1]$$

$$ \begin{align*} h[n] &= \alpha^n u[n]\\\\ H_{n,m} &= h[(n-m)_N] \end{align*} $$

%plot -s 560,200

N = 8;
nh = [0:N-1];
nx = [0:N-1];

x = zeros(size(nx))';
x(1) = 1;

h = zeros(size(nh))';
h = 0.8.^nh;

subplot(2,1,1), stem(nx,x,'filled'), ylabel('x[n]'), axis tight, ylim([0,1])
subplot(2,1,2), stem(nh,h,'filled'), ylabel('h[n]'), axis tight, ylim([0,1])

%plot -s 560,420

H = zeros(N,N);

for k = 0:N-1
    H(:,k+1) = circshift(h,k);
end

imagesc(H); colormap('jet'); axis('square'); axis off

3. Convolution

Convolution is defined as the integral of the product of the two functions after one is reversed and shifted

$$y[n] = \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m] = x[n]*h[n]$$

Output $y[n]$ came out by convolution of input $x[n]$ and system $h[n]$

• Commutative

$$x*h = h*x$$

$$ \begin{align*} y[n] &= \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m] =x[n]*h[n]\\ \\ &= \sum \limits_{k=-\infty}^{\infty} h[k]\,x[n - k] = \sum \limits_{k=-\infty}^{\infty} x[n - k]\,h[k]=h[n]*x[n]\\ \\ &\quad \;(k = n -m \, \implies \, m = n - k) \end{align*} $$

• signal = system

2.1. Convolution in MATLAB

$$\begin{align*} x[n]=\{3,11,7,&\;0,-1,4,2\}\quad \quad -3\leq n \leq 3\\ &\uparrow \\\\ h[n]=\{2,&\;3,0,-5,2,1\}\quad \quad -1\leq n \leq 4\\ &\uparrow \end{align*} $$

%plot -s 560,200

%% matlab command 'conv'
x = [3,11,7,0,-1,4,2]; 
h = [2,3,0,-5,2,1];    
y = conv(x,h);

stem(y,'filled'), ylim([-60 60])

% this is not correct

• You have to include conv_m function file in the path.

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

% Modified convolution routine for signal processing
% [y,ny] = conv_m(x,nx,h,nh)
%  y = convolution result
% ny = support of y
%  x = first signal on support nx
% nx = support of x
%  h = second signal on support nh
% nh = support of h

nyb = nx(1) + nh(1); 
nye = nx(length(x)) + nh(length(h));
ny = [nyb:nye];
y = conv(x,h);

$$\begin{align*} x[n]=\{3,11,7,&\;0,-1,4,2\}\quad \quad -3\leq n \leq 3\\ &\uparrow \\\\ h[n]=\{2,&\;3,0,-5,2,1\}\quad \quad -1\leq n \leq 4\\ &\uparrow \end{align*} $$

%plot -s 560,200

x = [3,11,7,0,-1,4,2];  nx = [-3:3];
h = [2,3,0,-5,2,1];     nh = [-1:4];

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

stem(ny,y,'filled'), axis tight, ylim([-60 60])

% this is what we want

%plot -s 560,600

subplot(3,1,1), stem(nx,x,'filled'); xlim([ny(1),ny(end)]);  ylabel('x')
subplot(3,1,2), stem(nh,h,'filled'); xlim([ny(1),ny(end)]);  ylabel('h')
subplot(3,1,3), stem(ny,y,'filled'); xlim([ny(1),ny(end)]);  ylabel('y')

2.2. Convolution For Infinite-Length Signals

(= Toeplitz Matrices)

$$\begin{align*}y[n] &= x[n]*h[n] = \sum \limits_{m=-\infty}^{\infty} h[n - m]\,x[m]\\\\ &= \cdots + \,h[n+2]\,x[-2] + h[n+1]\,x[-1] + h[n]\,x[0] + h[n-1]\,x[1] + h[n-2]\,x[2]\, + \cdots \end{align*}$$

It is an inner product of $h$ vectors and $x$

$$ \begin{align*} y[n] &= \begin{bmatrix} \cdots & \mid & \mid & \mid & \mid & \mid & \cdots\\ \cdots & h[n+2] & h[n+1] & h[n] & h[n-1] & h[n-2]\ & \cdots\\ \cdots & \mid & \mid & \mid & \mid & \mid & \cdots \end{bmatrix} \begin{bmatrix} \vdots\\ x[-2]\\ x[-1]\\ x[0]\\ x[1]\\ x[2]\\ \vdots \end{bmatrix}\\ \\ &= H x \end{align*} $$

Convolution is product of matrix $H$ and $x$

x = [3, 11, 7, 0, -1, 4, 2]';
h = [2, 3, 0, -5, 2, 1]';

hc = [h,zeros([length(x)-1,1])];
hr = [h(1),zeros([1,length(x)-1])];
H = toeplitz(hc,hr)

y = H*x

H =

     2     0     0     0     0     0     0
     3     2     0     0     0     0     0
     0     3     2     0     0     0     0
    -5     0     3     2     0     0     0
     2    -5     0     3     2     0     0
     1     2    -5     0     3     2     0
     0     1     2    -5     0     3     2
     0     0     1     2    -5     0     3
     0     0     0     1     2    -5     0
     0     0     0     0     1     2    -5
     0     0     0     0     0     1     2
     0     0     0     0     0     0     1


y =

     6
    31
    47
     6
   -51
    -5
    41
    18
   -22
    -3
     8
     2

2.3. Convolution For Finite-Length Signals

Circular Convolution

$$y[n] = x[n] \otimes h[n] = \sum_{m = 0}^{N-1} h \left[(n-m)_N \right]\,x[m]$$

• Circularly time reverse the impulse response $h$ and circularly shift it $n$ time steps to the right (delay)

• Compute the inner product between the shifted impulse response and the input vector $x$

%plot -s 560,200

N = 8;
n = 0:N-1;

h = [0 1 2 3 4 5 6 7]';

stem(n,h,'filled');

%% Create a circulent matrix H based on h

H = zeros(N,N);
for k = 0:N-1
    H(:,k+1) = circshift(h,k);
end
H

H =

     0     7     6     5     4     3     2     1
     1     0     7     6     5     4     3     2
     2     1     0     7     6     5     4     3
     3     2     1     0     7     6     5     4
     4     3     2     1     0     7     6     5
     5     4     3     2     1     0     7     6
     6     5     4     3     2     1     0     7
     7     6     5     4     3     2     1     0

%plot -s 800,300

subplot(1,2,1), imagesc(H); colormap('jet'); axis('square'); axis off
subplot(1,2,2), imagesc(h); colormap('jet'); axis equal; axis off

x = [1,1,0,0,0,0,0,0]';
y1 = H*x

y1 =

     7
     1
     3
     5
     7
     9
    11
    13

x = [1,1,0,0,0,0,0,0];
h = [0,1,2,3,4,5,6,7];

y3 = cconv(x,h,8)

y3 =

    7.0000    1.0000    3.0000    5.0000    7.0000    9.0000   11.0000   13.0000

x1 = [1,2,2]; 
x2 = [1,2,3,4];

y = cconv(x1,x2,4)

% 5-point circular convolution
y = cconv(x1,x2,5)

% 6-point circular convolution
y = cconv(x1,x2,6)

y =

    15    12     9    14


y =

     9     4     9    14    14


y =

     1     4     9    14    14     8

%plot -s 560,500

N = 8;
xn = 0:N-1;
x = [1 0 -1 0 1 0 -1 0];
h = [0 1 2 3 4 0 0 0];
y = cconv(x,h,N);

subplot(3,1,1), stem(xn,x,'filled'), ylim([-2 4]), ylabel('x')
subplot(3,1,2), stem(xn,h,'filled'), ylim([-2 4]), ylabel('h')
subplot(3,1,3), stem(xn,y,'filled'), ylim([-2 4]), ylabel('y')

%plot -s 560,600

N = 8;
theta2 = (0:1/N:1)*2*pi;
theta = theta2(1:end-1);

re = cos(theta);
im = sin(theta);

subplot(3,1,1),  stem3(re,im,x,'filled');   hold on;
plot(re,im,'k')
quiver(0,0,1,0,'k','maxheadsize',0.5)
quiver(0,0,0,1,'k','maxheadsize',0.5), hold off
xlabel('Real','fontsize',8), ylabel('Imaginary','fontsize',8), zlabel('x[n]')
zlim([-2,4])

subplot(3,1,2),  stem3(re,im,h,'filled');   hold on;
plot(re,im,'k')
quiver(0,0,1,0,'k','maxheadsize',0.5)
quiver(0,0,0,1,'k','maxheadsize',0.5), hold off
xlabel('Real','fontsize',8), ylabel('Imaginary','fontsize',8), zlabel('h[n]')
zlim([-2,4])

subplot(3,1,3),  stem3(re,im,y,'filled');   hold on;
plot(re,im,'k')
quiver(0,0,1,0,'k','maxheadsize',0.5)
quiver(0,0,0,1,'k','maxheadsize',0.5), hold off
xlabel('Real','fontsize',8), ylabel('Imaginary','fontsize',8), zlabel('y[n]')
zlim([-2,4])

4. Convolution Examples

4.1. Denoising

%plot -s 560,800

N = 50;
n = 0:N-1;
s = sin(pi/N*n) .* [ones(1,N/2), -ones(1,N/2)];

subplot(4,1,1)
stem(s,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('piecewise smooth signal','fontsize',8)

% add noise to the signal
x = s + 0.1*randn(1,N);

subplot(4,1,2)
stem(x,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('piecewise smooth signal + noise','fontsize',8)

% construct moving average filter impulse response of length M
M = 5;
h = ones(1,M)/M;

subplot(4,1,3)
stem(h,'b','Marker','none','LineWidth',2);
axis([1 N 0 1])
title('impulse response','fontsize',8)

% convolve noisy signal with impulse response
y = cconv(x,h,N);

subplot(4,1,4)
stem(y,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('output','fontsize',8)

4.2. Edge Detection

%plot -s 560,500

N = 50;
n = 0:N-1;
s = sin(pi/N*n) .* [ones(1,N/2), -ones(1,N/2)];

% add noise to the signal
x = s + 0.1*randn(1,N);

subplot(311)
stem(x,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('piecewise smooth signal + noise','fontsize',8)

% haar wavelet edge detector
h = (1/4)*[-ones(1,2), ones(1,2)];

subplot(312)
stem(h,'b','Marker','none','LineWidth',2);
axis([1 N -0.5 0.5])
title('Haar wavelet edge detector','fontsize',8)

% % convolve noisy signal with impulse response
y = cconv(x,h,N);

subplot(313)
stem(y,'b','Marker','none','LineWidth',2);
axis([1 N -1.5 1.5])
title('output','fontsize',8)

4.3. Convolution on Audio

• ThinkDSP by Allen Downey at Olin (http://nbviewer.ipython.org/github/AllenDowney/ThinkDSP/blob/master/code/scipy2015_demo.ipynb)

• Download wave file 1

• Download wave file 2

[x, Fs] = audioread([pwd,'\data_files\violin_origional.wav']);

x = x/max(x);  % normalize
sound(x, Fs);  % play a wave file with sampling rate Fs

%plot -s 560,200

% plot wave file
xn = (1:length(x))/Fs;
plot(xn,x), xlabel('time in sec')

% impulse response in a closed room (by gunshot)

[h, Fs] = audioread([pwd,'\data_files\gunshot.wav']);
h = h(:,1); % stereo to mono

sound(h, Fs);
plot((1:length(h))/Fs, h), xlabel('sec')

% image how the music played in a closed room sounds like

y = cconv(x,h);
y = y/max(y);

sound(y, Fs);

%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')