Make Your Own Convolution

by

H. Joel Trussell
Dept. of Electrical and Computer Engineering
North Carolina State University
Raleigh, NC  27695-7911

e-mail:  hjt@ecehjt.ncsu.edu

This notebook let you create your own function and convolution kernel with which to demonstrate convolution.

Initialization

makeFunction[x_] := Block[{ m,n,y},

	m = Length[x];
	n = 1;
	y = {x[[1]]};

	For[i=2,i<=m,i++,

		If[x[[i,1]] - y[[n,1]] > 0.01,

			(* then *) y = Append[y,x[[i]]];
			n = n+1;

		];

	];
	Return[y];
];

Make Your Function

Enter your function on the graph below. Click once in the graph area to set the frame around the axes. To draw you function, hold down the Command key and the left button on the mouse as you move the mouse. The points will appear on the graph. You can repeat the process until you have the graph you want. Next, press <Command>-c (hold the Command key while pressing the c key) to copy the point pairs that define your function. A function f(x) will be defined using those pairs. This calls for you to alter the command in the Mathematica file by deleting any previous lists of points after inserting the
list you have just copied. Insert the list in the place indicated below.

ListPlot[{{-10,-5},{10,5}}]

(* insert list of point pairs after the = sign for the *)
(* variable x below *)
x = {{-2.954321, 0.045341}, {-2.055959, 0.590436}, 
{-1.23246, -0.015225}, {-0.858143, -0.74202}, 
{-0.109508, -1.771645}, {0.189947, -1.953344}, 
{0.788855, -2.195609}, {1.3129, -1.468814}, 
{1.836945, -0.74202}, {2.435853, 0.105907}, 
{3.558805, 0.893268}, {3.85826, 0.832702}, 
{4.157714, 0.529871}, {4.532032, 0.408738}, 
{4.681759, 0.105907}};
fs  = makeFunction[x]; (*eliminate points with duplicate *)
						(* abscissa values  *)
f = Interpolation[fs,InterpolationOrder->1]; (* defines the function f(x) *)
fxmin = Min[Table[fs[[i,1]],{i,Length[fs]}]];
fxmax = Max[Table[fs[[i,1]],{i,Length[fs]}]];
ListPlot[Table[{t,f[t]},{t,fxmin,fxmax,(fxmax-fxmin)/128}],
	PlotJoined->True,
	PlotLabel -> "Function in Time Domain"
]

Make Your Kernel

Just as you defined your function, now define the kernel or impulse response that characterizes your system. Notice that the scale has been changed. This reflects the fact that most impulse responses are much shorter that the signals that are processed.

ListPlot[{{-2,-1},{2,1}}]

(* insert list of point pairs after the = sign for the *)
(* variable x below *)
x = {{0.014924, 0.021111}, {0.029847, 0.081478}, 
{0.029847, 0.129771}, {0.029847, 0.190138}, 
{0.029847, 0.250504}, {0.029847, 0.298797}, 
{0.029847, 0.359164}, {0.029847, 0.407457}, 
{0.029847, 0.467823}, {0.029847, 0.516116}, 
{0.029847, 0.588556}, {0.059694, 0.660996}, 
{0.059694, 0.709289}, {0.08954, 0.757582}, 
{0.134311, 0.709289}, {0.164158, 0.660996}, 
{0.194005, 0.612702}, {0.223851, 0.564409}, 
{0.253698, 0.516116}, {0.298468, 0.467823}, 
{0.343239, 0.41953}, {0.402932, 0.371237}, 
{0.47755, 0.31087}, {0.537243, 0.27465}, 
{0.596937, 0.238431}, {0.656631, 0.202211}, 
{0.761095, 0.165991}, {0.820788, 0.117698}, 
{0.895405, 0.093551}, {0.970022, 0.057331}, 
{1.059563, 0.033185}, {1.119256, 0.009038}, 
{1.17895, -0.003035}, {1.238644, -0.003035}};
hs  = makeFunction[x]; (*eliminate points with duplicate *)
						(* abscissa values  *)
h = Interpolation[hs,InterpolationOrder->1]; (* defines the function f(x) *)
hxmin = Min[Table[hs[[i,1]],{i,Length[hs]}]];
hxmax = Max[Table[hs[[i,1]],{i,Length[hs]}]];
ListPlot[Table[{t,h[t]},
	{t,hxmin+0.001,hxmax-0.001,(hxmax-hxmin-0.002)/128}],
	PlotJoined->True,
	PlotLabel -> "Impulse Response in Time Domain"
]

Compute Your Convolution

The following commands set up the product for the convolutional integral and perform the  integration. For the Interpolation function to yield the proper values in the range of the integral, the function must be defined over the entire range of the integral. We define it zero in the regions where you have not defined it.

If [hxmin < 0,
	fsx = Join[{{fxmin+hxmin,0},{(2 fxmin+hxmin)/2,0}},fs],
	(* else *) fsx = fs];
If [hxmax > 0,
	fsx = Join[fsx, {{(2 fxmax + hxmax)/2,0},{fxmax+hxmax,0}}]];
f = Interpolation[fsx,InterpolationOrder->1];
If [fxmin < 0,
	hsx = Join[{{fxmin+hxmin,0},{(fxmin+(2 hxmin))/2,0}},hs],
	(* else *) hsx = hs];
If [fxmax > 0,
	hsx = Join[hsx, {{(fxmax + (2hxmax))/2,0},{fxmax+hxmax,0}}]];
h = Interpolation[hsx,InterpolationOrder->1];

lowerLim = fxmin+hxmin;
upperLim = fxmax+hxmax;
deltaT = (upperLim-lowerLim)/255;
df = Table[f[t],{t,lowerLim,upperLim,deltaT}];
dh = Table[h[t],{t,lowerLim,upperLim,deltaT}];
If [lowerLim < 0, nshift = Floor[-lowerLim/deltaT], nshift = 0];
dhs = dh;
For[i=0,i<Length[dh],i++,dhs[[i+1]] = dh[[Mod[i+nshift,Length[dh]]+1]]];
CircularFFT[x_,h_] := Block[{Nx,Nh,Ny,xpad,hpad,
						FFTxpad,FFThpad,FFTy},
	Nx = Length[x];
	Nh = Length[h];
	If[Nx != Nh, Print["Signals are not same length"];
			Return];
	Ny = Nx;
	FFTx = Fourier[x];
	FFTh = Fourier[h];
	FFTy = FFTx FFTh;
	y = N[Sqrt[Ny]] InverseFourier[FFTy];
	Return[Re[y]];
];
dy = CircularFFT[df,dhs];

Plot Results of convolution

timeAxis = Table[t,{t,lowerLim,upperLim,deltaT}];
Show[ GraphicsArray[{ {
  Plotfunction = ListPlot[
  	Transpose[{timeAxis,df}],
    PlotJoined->True,
	PlotLabel -> "Function in Time Domain",
	DisplayFunction->Identity,
	PlotRange->All],
  Plotkernel = ListPlot[
  	Transpose[{timeAxis,dh}],
    PlotJoined->True,
	PlotLabel -> "Impulse Response in Time Domain",
	DisplayFunction->Identity,
	PlotRange->All]},{
  Plotresult = ListPlot[
  	Transpose[{timeAxis,dy}],
    PlotJoined->True,
	PlotLabel -> "Output of Convolution",
	DisplayFunction->Identity],
  NullPlot = ListPlot[{{-10,-5},{10,5}},
    DisplayFunction->Identity]}}],
  	DisplayFunction :> $DisplayFunction
]