// convolution.sq - convolution and correlation of signals
// Version 2.5, 1.8.2003
// Copyright 1998-2003, Calerga.
//
//   Disclaimer
//
// Absolutely no guarantee is made about the content of this file, the
// underlying algorithms, their suitability with respect to any particular
// application, or the functionning of the code with any software. The user
// supports the whole responsibility of the reading and the use of this file
// and the use of the results which might be obtained with it.
//
// Permission is granted to licensees of SQ to use, modify, and distribute
// modified versions of this file or new files which borrow some code from
// it to other licensees of SQ, provided that its original author (Yves Piguet)
// is acknowledged and that it is made clear that the new file is a derivated
// work. The licensee assumes all the legal consequences which might result
// from the use and the distribution of the derivative work.
//
// The "data" block at the end of this file, if there is one, was added
// after the release of this file and is even more outside of the scope of
// the author's responsibility.

title "Convolution product"

version
{@
Version 2.5, 1st August 2003

Copyright 1998-2003, Calerga Sarl.
@}

help
{@
This SQ file illustrates how convolution and correlation products are
calculated.

The convolution product of two signals x(t) and y(t), denoted x(t)*y(t),
is defined as

x(t) * y(t) = integral(tau = -inf..inf, x(tau) y(t-tau) dtau)

At the beginning, the convolution product of two exponential
signals is performed. The upper figure shows the integral for a given value
of t, represented as a red circle in the result displayed in the lower
figure. The integral, shown as a yellow region, is the product of x and a
flipped y shifted by t. You can change t with the mouse either by sliding
the whole flipped y in the upper figure, or by moving the red circle on the
result. You can change the signals in the Settings menu.
@}

variable fn		// 0 = squares, 1 = exps, 10 = autocorrelation of exp
variable t		// shift of y

init (fn, t) = init

menu "Squares" _checkmark(fn == 0) fn = 0
menu "Exponentials" _checkmark(fn == 1) fn = 1
menu "Delayed Exponential" _checkmark(fn == 2) fn = 2
menu "Echo" _checkmark(fn == 3) fn = 3
menu "Step Response" _checkmark(fn == 4) fn = 4
separator
menu "Autocorrelation of Exponentials" _checkmark(fn == 10) fn = 10
menu "Signal Detection" _checkmark(fn == 11) fn = 11

figure "First Operand"
	draw drawX(fn)

figure "Second Operand"
	draw drawY(fn)

figure "Operands"
	draw drawOp(fn, t)
	mousedrag (t, _msg) = dragOp(t, _x, _x1)
	mouseover _msg = overOp(_id, fn, t)

figure "Result"
	draw drawResult(fn, t)
	mousedrag (t, _msg) = dragOp(t, _x0, _x1)
	mouseover (_msg, _cursor) = overResult(_id, _x0)

functions
{@
function (fn, t) = init
	
	fn = 1;
	t = 1;
	subplots('Operands\nResult');

function drawX(fn)
	
	scale([-5, 5, 0, 1.2]);
	sc = scale;
	tmin = sc(1);
	tmax = sc(2);
	switch fn
		case 0
			(t, x) = sigSquare(tmin, tmax, 0, 1);
			plot(t, x);
		case 1
			(t, x) = sigExp(tmin, tmax, 0, 2, false);
			plot(t, x);
		case 2
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x);
		case 3
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x);
		case 4
			(t, x) = sigStep(tmin, tmax, 0, false);
			plot(t, x);
		case 10
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x);
		case 11
			randn('s', 0);
			t = tmin + (tmax - tmin) * (0:300) / 300;
			x = (t >= 0 & t <= 1) + 0.5 * randn(size(t));
			scale([-5, 5, -1, 2]);
			plot(t, x);
	end

function drawY(fn)
	
	scale([-5, 5, 0, 1.2]);
	sc = scale;
	tmin = sc(1);
	tmax = sc(2);
	switch fn
		case 0
			(t, x) = sigSquare(tmin, tmax, 0, 1);
			plot(t, x);
		case 1
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x);
		case 2
			plot([1, 1], [0, sc(4)]);						% y = dirac(t-1)
		case 3
			plot([0, 0], [0, sc(4)]);						% y = dirac(t) + dirac(t-1)
			plot([1, 1], [0, sc(4)]);
		case 4
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x);
		case 10
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x);
		case 11
			(t, x) = sigSquare(tmin, tmax, 0, 1);
			plot(t, x);
	end

function drawOp(fn, tau)
	
	scale([-5, 5, 0, 1.2]);
	sc = scale;
	tmin = sc(1);
	tmax = sc(2);
	t = tmin + (tmax - tmin) * (0:1000) / 1000;
	line([1, 0], tau, 'c', 100);
	lgnd = 'f1(t)\nf2(tau-t)\n(f1*f2)(tau)';
	switch fn
		case 0
			plot([tmin,t,tmax], [0,(t>=0&t<=1).*(tau-t>=0&tau-t<=1),0], 'yf');
			(t, x) = sigSquare(tmin, tmax, 0, 1);
			plot(t, x, '', 1);
			(t, x) = sigSquare(tmin, tmax, tau - 1, tau);
			plot(t, x, 'r', 2);
		case 1
			plot([tmin,t,tmax], [0,exp(-t/2).*(t>=0).*exp(t-tau).*(tau-t>=0),0], 'yf');
			(t, x) = sigExp(tmin, tmax, 0, 2, false);
			plot(t, x, '', 1);
			(t, x) = sigExp(tmin, tmax, tau, 1, true);
			plot(t, x, 'r', 2);
		case 2
			plot((tau-1)+[-0.05,0.05,0.05,-0.05], exp(-tau+1)*(tau-1>=0)*[0,0,1,1], 'yf');
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x, '', 1);
			plot([tau, tau]-1, [0, sc(4)], 'r', 2);			% y = dirac(t-1)
		case 3
			plot(tau+[-0.05,0.05,0.05,-0.05], exp(-tau)*(tau>=0)*[0,0,1,1], 'yf');
			plot((tau-1)+[-0.05,0.05,0.05,-0.05], exp(-tau+1)*(tau-1>=0)*[0,0,1,1], 'yf');
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x, '', 1);
			plot([tau, tau], [0, sc(4)], 'r', 2);			% y = dirac(t) + dirac(t-1)
			plot([tau, tau]-1, [0, sc(4)], 'r', 2);
		case 4
			plot([tmin,t,tmax], [0,exp(t-tau).*(t>=0&t<=tau),0], 'yf');
			(t, x) = sigStep(tmin, tmax, 0, false);
			plot(t, x, '', 1);
			(t, x) = sigExp(tmin, tmax, tau, 1, true);
			plot(t, x, 'r', 2);
		case 10
			plot([tmin,t,tmax], [0,exp(-t).*(t>=0).*exp(-t+tau).*(t-tau>=0),0], 'yf');
			(t, x) = sigExp(tmin, tmax, 0, 1, false);
			plot(t, x, '', 1);
			(t, x) = sigExp(tmin, tmax, tau, 1, false);
			plot(t, x, 'r', 2);
			lgnd = 'f(t)\nf(t+tau)\nR_ff(tau)';
		case 11
			scale([-5, 5, -1, 2]);
			randn('s', 0);
			t = tmin + (tmax - tmin) * (0:300) / 300;
			x = (t >= 0 & t <= 1) + 0.5 * randn(size(t));
			plot([tmin, t, tmax], [0, x.*(t >= tau & t <= tau + 1), 0], 'yf');
			plot(t, x, '', 1);
			(t, x) = sigSquare(tmin, tmax, tau, tau + 1);
			plot(t, x, 'r', 2);
			lgnd = 'f1(t)\nf2(t+tau)\nR_f1f2(tau)';
	end
	
	legend(lgnd, 'k_r_yf');

function drawResult(fn, tau)
	
	scale([-5, 5]);
	sc = scale;
	tmin = sc(1);
	tmax = sc(2);
	t = tmin + (tmax - tmin) * (0:1000) / 1000;
	lgnd = '(f1*f2)(t)\n(f1*f2)(tau)';
	switch fn
		case 0
			plot(t, (t>=0).*t-(t>=1).*(2.*t-2)+(t>=2).*(t-2));
			plot(tau, (tau>=0).*tau-(tau>=1).*(2.*tau-2)+(tau>=2).*(tau-2), 'or', 1);
		case 1
			plot(t, convExp(1, 2, t));
			plot(tau, convExp(1, 2, tau), 'or', 1);
		case 2
			plot(t, exp(-t+1) .* (t-1 >= 0));
			plot(tau, exp(-tau+1) .* (tau-1 >= 0), 'or', 1);
		case 3
			plot(t, exp(-t) .* (t >= 0) + exp(-t+1) .* (t-1 >= 0));
			plot(tau, exp(-tau) .* (tau >= 0) + exp(-tau+1) .* (tau-1 >= 0), 'or', 1);
		case 4
			plot(t, (1 - exp(-t)) .* (t >= 0));
			plot(tau, (1 - exp(-tau)) .* (tau >= 0), 'or', 1);
		case 10
			plot(t, exp(-3 * abs(t)) / 2);
			plot(tau, exp(-3 * abs(tau)) / 2, 'or', 1);
			lgnd = 'R_ff(t)\nR_ff(tau)';
		case 11
			randn('s', 0);
			t = tmin + (tmax - tmin) * (0:ceil(300 * (tmax + 1 - tmin) / (tmax - tmin))) / 300;
			x = (t >= 0 & t <= 1) + 0.5 * randn(size(t));
			delta = round(300 / (tmax - tmin));
			cor = (cumsum(x((1+delta:length(t)))) - cumsum(x(1:length(t)-delta))) / delta;
			plot(t(1:length(cor)), cor);
			(dmy, ix) = min(abs(tau - t));
			if ix <= length(cor)
				plot(tau, cor(ix), 'or', 1);
			end
			lgnd = 'R_f1f2(t)\nR_f1f2(tau)';
	end
	legend(lgnd, '_kor');

function (t, msg) = dragOp(t, x, x1)
	
	if isempty(x)
		cancel;
	end
	t = t + x1 - x;
	msg = sprintf('t = %g', t);

function msg = overOp(id, fn, t)
	
	if isempty(id)
		cancel;
	end
	switch id
		case 1
			switch fn
				case 0
					msg = 'x(t) = 1 if 0 <= t <= 1, 0 otherwise';
				case 1
					msg = 'x(t) = exp(-t/2) for t >= 0';
				case 2
					msg = 'x(t) = exp(-t) for t >= 0';
				case 3
					msg = 'x(t) = exp(-t) for t >= 0';
				case 4
					msg = 'x(t) = exp(-t) for t >= 0';
				case 10
					msg = 'x(t) = exp(-t) for t >= 0';
				case 11
					msg = 'x(t) = (1 if 0 <= t <= 1, 0 otherwise) + noise';
			end
		case 2
			switch fn
				case 0
					msg = 'y(t) = 1 if 0 <= t <= 1, 0 otherwise';
				case 1
					msg = 'y(t) = exp(-t) for t >= 0';
				case 2
					msg = 'y(t) = dirac(t-1)';
				case 3
					msg = 'y(t) = dirac(t) + dirac(t-1)';
				case 4
					msg = 'y(t) = step(t)';
				case 10
					msg = 'y(t) = exp(-t) for t >= 0';
				case 11
					msg = 'y(t) = 1 if 0 <= t <= 1, 0 otherwise';
			end
		case 100
			msg = sprintf('t = %g', t);
	end

function (msg, cursor) = overResult(id, x0)
	
	if id ~== 1
		cancel;
	end
	msg = sprintf('t = %g', x0);
	cursor = true;

function (t, x) = sigSquare(tmin, tmax, stepmin, stepmax)
	
	if stepmin > tmax | stepmax < tmin
		t = [tmin, tmax];
		x = [0, 0];
	elseif stepmin < tmin && stepmax > tmax
		t = [tmin, tmax];
		x = [1, 1];
	elseif stepmin < tmin
		t = [tmin, stepmax, stepmax, tmax];
		x = [1, 1, 0, 0];
	elseif stepmax > tmax
		t = [tmin, stepmin, stepmin, tmax];
		x = [0, 0, 1, 1];
	else
		t = [tmin, stepmin, stepmin, stepmax, stepmax, tmax];
		x = [0, 0, 1, 1, 0, 0];
	end

function (t, x) = sigStep(tmin, tmax, t0, backward)
	
	if t0 >= tmin && t0 <= tmax
		t = [tmin, t0, t0, tmax];
		if backward
			x = [1, 1, 0, 0];
		else
			x = [0, 0, 1, 1];
		end
	else
		t = [tmin, tmax];
		x = [0, 0] + (backward && t0 > tmax | ~backward && t0 < tmin);
	end

function (t, x) = sigExp(tmin, tmax, t0, tau, backward)
	
	if backward
		if t0 > tmax
			t = tmin + (tmax - tmin) * (0:100) / 100;
			x = exp((t - t0) / tau);
		elseif t0 < tmin
			t = [tmin, tmax];
			x = [0, 0];
		else
			t = [tmin+(t0-tmin)*(0:100)/100, t0, tmax];
			x = [exp((tmin+(t0-tmin)*(0:100)/100-t0)/tau), 0, 0];
		end
	else
		if t0 < tmin
			t = tmin + (tmax - tmin) * (0:100) / 100;
			x = exp((t0 - t) / tau);
		elseif t0 > tmax
			t = [tmin, tmax];
			x = [0, 0];
		else
			t = [tmin, t0, t0+(tmax-t0)*(0:100)/100];
			x = [0, 0, exp((t0-(t0+(tmax-t0)*(0:100)/100))/tau)];
		end
	end

function y = convExp(tau1, tau2, t)
	
	if tau1 == tau2
		y = t .* exp(-t ./ tau1);
	else
		y = tau1 * tau2 / (tau1 - tau2) .* (exp(-t / tau1) - exp(-t / tau2)) .* (t >= 0);
	end

@}