%?ving 7  Maximum likelihood estimator av bevegelse

%% Modell for autokorrelasjon av x , bilde til bilde

f0=2e6;
c=1540;
lambda=c/f0;
As=1; An=0.3;

dzTrue=0.1e-3;

%simulering av signal x1 x2

s1=As*(randn+i*randn);%signal fra bilde 1  Kompleks Gaussisk
s2=s1*exp(i*4*pi*dzTrue/lambda);%signal fra bilde 2.  Fase-skift p? grunn av radiell forskyvning dzTrue

x1=s1+An*(randn+i*randn);%Additiv Kompleks Gaussisk  st?y
x2=s2+An*(randn+i*randn);%Additiv Kompleks Gaussisk  st?y.  Ukorrelert fra bilde til bilde
%return;
maxAmp=max(abs(x1),abs(x2));
maxAmp=2;
plot([0,x1],'r');hold on;plot([0,x2],'b');
axis equal;grid;axis(maxAmp*[-1,1,-1,1]);
hold off;
%return;

%% Numerisk beregning av ML estimat:

dz=linspace(-lambda/4,lambda/4,200);%Kandidater for radiell forskyvning

x=[x1,x2];%signal vektor

for k=1:length(dz),
    R0=As^2 + An^2;% 0.5*< x1'*x1 + x2'*x2 >   (  < > betyr forventningsverdi)
    R1=As^2*exp(i*4*pi*dz(k)/lambda);% < x1'*x2 >
    C=[R0,R1;conj(R1),R0];%Kovarians matrise
    p=1/sqrt(det(C))*exp(-0.5*x*inv(C)*x');% sans. tetthet for x, gitt dz(k) forskyvning
    l(k)=p; % likelihood funksjon
end;

plot(dz,l);grid;
[val,ind]=max(l);%finner indeks for max verdi av l
dzEst=dz(ind)% Og tilsvarende dz verdi, som er ML estimatet av dz


