% M03-de+corr.diary
% Running an LCCDE, plus cross-correlation
% 2007-09-18

% Evaluating an LCCDE algorithmically
% y[n] + y[n-1] - 6y[n-2] = x[n]
% with y[-2] = -1, y[-1] = 1, x[n] = 8.u[n]

% Calculate for n = 0..4, but include -2 and -1
nn = -2:4;

% So setup initial y values; zero is index 3
zz = 3;
y(zz+ -2) = -1;
y(zz+ -1) = 1;
% Now evaluate the LCCDE step by step
x = zeros(1,7);
x(zz+ [0:4]) = 8;

for n = 0:4; y(zz+ n) = x(zz+ n) - y(zz+ n-1) + 6*y(zz+ n-2); end
y
% y =
%     -1     1     1    13     1    85   -71

% Compare to algebraic solution
l1 = -3;
l2 = 2;
beta = -2;
a1 = -1.8;
a2 = 4.8;
yt = a1*(l1.^nn) + a2*(l2.^nn) + beta
% yt =
%    -1.0000    1.0000    1.0000   13.0000    1.0000   85.0000  -71.0000
nn
% nn =
%     -2    -1     0     1     2     3     4

% In this case, yt[-2, -1] *do* match ICs, but not always true

% Use built-in filter() command to calculate impulse response
% (you can't (easily) use filter() if there are initial conditions)
% x DE coefficients {p_k} = {1}, y DE coefficients {d_k} = {1,1,-6}
ir = filter([1], [1 1 -6], [1 0 0 0 0 0 0 0])
% ir =
%   Columns 1 through 4 
%            1          -1           7         -13
%   Columns 5 through 8 
%           55        -133         463       -1261

% Compare to closed form
h = 0.6*((-3).^[0:7]) + 0.4*(2.^[0:7])
% h =
%    1.0e+03 *
%   Columns 1 through 5 
%     0.0010   -0.0010    0.0070   -0.0130    0.0550
%   Columns 6 through 8 
%    -0.1330    0.4630   -1.2610

% different formatting, but the same values...


% Cross-correlation visualization, to find a pulse in noise
% Basic pulse
p = [-1 3 2 1];
% Noise sequence
n = randn(1,20);
% Make the signal be the noise plus pulse at time 10
x = n;
x(10+[1:4]) = x(10+[1:4]) + p;
% Look at them
subplot(311)
stem(10:13,p);
axis([-2 20 -5 5]); grid on;
title('Pulse p');
subplot(312)
stem([0:19],n);
axis([-2 20 -5 5]); grid on;
title('Noise n');
subplot(313)
stem([0:19],x);
axis([-2 20 -5 5]); grid on;
title('x = noise + pulse');

% Calculate the cross-correlation step by step; max should be at 10
corrviz(p,x)