% eda07_06 Script supporting CribSheet 7.2
% Least squares filter estimation of filter m
% obeying g*m=d
clear all;
global g H;

% definition of time series
N=101;
Dt = 1;
t = Dt*[0:N-1]';
tmin = t(1);
tmax = t(end);

% make some synthetic data g
g = random('exp',1,N,1);

% true filter is decaying exponential
M=10;
c=5*Dt;
mtrue=0.3*exp(-t(1:M)/c);
dtrue = conv(g,mtrue);
dtrue = dtrue(1:N);

% observed data is true data plus noise
sigmad = 0.1;
d = dtrue + random('Normal',0,sigmad,N,1);

% plot g, d
figure(1);
clf;
subplot(2,1,1)
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, -5, 5] );
plot( t, g, 'k-', 'LineWidth', 2 );
xlabel('t');
ylabel('g');
subplot(2,1,2)
set(gca,'LineWidth',2);
hold on;
axis( [tmin, tmax, -5, 5] );
plot( t, d, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('d');

% plot filter m
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [tmin, t(M), 0, 1] );
plot( t(1:M), mtrue, 'k-', 'LineWidth', 3 );
xlabel('t');
ylabel('m');

% prior infoirmation of smoothness
epsilon = 0.001;
e = epsilon*ones(M,1)/(Dt^2);
H = spdiags([e,-2*e,e],[-1:1],M-2,M);
h = zeros(M-2,1);
K=M-2;

% rhs of least squares equations
L=N+K;
f(1:N)=d;
f(N+1:N+K)=h;
temp=conv(flipud(g),d);
FTfa = temp(N:N+M-1);
FTfb = H'*h;
FTf=FTfa+FTfb;

% least squares solution
mest=bicg(@filterfun,FTf,1e-10,3*L);

% predicted data
dpre = conv(g,mest);
dpre = dpre(1:N);
e = d-dpre;
E=e'*e;

% plot predicted data
figure(1);
subplot(2,1,2)
set(gca,'LineWidth',2);
hold on;
plot( t, dpre, '-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8] );

% plot estimated filtes
figure(2);
set(gca,'LineWidth',2);
hold on;
axis( [tmin, t(M), 0, 1] );
plot( t(1:M), mest, '-', 'LineWidth', 2, 'Color', [0.8,0.8,0.8] );

