% compute impulse response with least squares
clear all;

% definitionof time with N points
N = 201;
tmin = 0;
tmax = 50;
dt = (tmax-tmin)/(N-1);
tn = tmin+dt*[0:N-1]';

% definition of time with fewer points
M=31;
tm = tmin+dt*[0:M-1]';

% the true impulse response is a decaying exponential
c=1;
gn = exp(-c*tn);

% matrix form of the impulse response
G=zeros(N,N);
for p = [1:N]
    G(p:N,p) = gn(1:N-p+1);
end
G=dt*G;

% the heating is the product of a sine and an exponential, with additive
% constant
f=0.2;
d=0.1;
hn = 2+sin(2*pi*f*tn).*exp(-d*tn);

% temperature is impulse response acting on heating
zn = G*hn;

% plot impulse response, heating, temperature
figure(1);
clf;
subplot(3,1,1);
plot(tn,gn);
xlabel('time, t');
ylabel('impulse response');
hold on;
subplot(3,1,2);
plot(tn,hn);
xlabel('time, t');
ylabel('heating');
hold on;
subplot(3,1,3);
plot(tn,zn);
xlabel('time, t');
ylabel('temperature');
hold on;

% add random noise to heating and plot
figure(2);
clf;
h2n = hn + 0.2*(rand(N,1)-0.5);
subplot(2,1,1);
plot(tn,h2n,'r');
hold on;
plot(tn,hn,'k');
xlabel('time, t');
ylabel('heating+noise');
hold on;

% add random noise to temperature and plot
z2n = zn + 0.5*(rand(N,1)-0.5);
subplot(2,1,2);
plot(tn,z2n,'r');
hold on;
plot(tn,zn,'k');
xlabel('time, t');
ylabel('temperature+noise');
hold on;

% matrix form of noisy heating
H=zeros(N,N);
for p = [1:N]
    H(p:N,p) = h2n(1:N-p+1);
end
H=dt*H;

% estimate impulse response with matrix inversion (yuck) and plot
gnest = inv(H)*z2n;
figure(3);
clf;
plot(tn,gn,'r');
hold on;
plot(tn,gnest,'k');

% matrix form of noisy heating, but for smaller number of model parameters
Hm=zeros(N,M);
for p = [1:M]
    Hm(p:N,p) = h2n(1:N-p+1);
end
Hm=dt*Hm;

% estimate impulse response by least squares
figure(4);
clf;
epsilon2=10.0;
gmest = inv(Hm'*Hm+epsilon2*dt*dt*eye(M))*Hm'*z2n;
plot(tn,gn,'r');
hold on;
plot(tm,gmest,'k');

% coefficients for second-derivative damping
a=1.0;
b=(-2.0);
c=1.0;

% set up D matrix
D=zeros(M-2,M);
for p = [1:M-2]
    D(p,p) = a;
    D(p,p+1) = b;
    D(p,p+2) = c;
end
D=dt*D;

% generalized least squar
figure(5);
clf;
epsilon2=100.0;
gmest2 = inv(Hm'*Hm+epsilon2*D'*D)*Hm'*z2n;
plot(tn,gn,'r');
hold on;
plot(tm,gmest2,'k');