
Ta=-5;
Tb=5;
N=10001;
t=linspace(Ta,Tb,N);
ns=(N-1)/(Tb-Ta);
Ts=1;
tau1=0.5*rand;
tau2=0.5*rand;
t1=round(tau1*ns);
t2=round(tau2*ns);
eta1=0.8+0.2*rand;
eta2=0.8+0.2*rand;
phi1=-pi/4+pi/2*rand;
phi2=-pi/4+pi/2*rand;
c1=round(rand);
c2=round(rand);

if c1==1
    d1=1;
elseif c1==0
    d1=-1;
end
if c2==1
    d2=1;
elseif c2==0
    d2=-1;
end

sat=zeros(1,N);
sat(4501+t1:6502+t1)=1;

sbt=zeros(1,N);
sbt(4501+t2:5501+t2)=1;
sbt(5502+t2:6502+t2)=-1;
phi=2*pi*randn(1,N);
wt=0.1*exp(i*phi).*randn(1,N);

xt=eta1*exp(i*phi1)*d1*sat+eta2*exp(i*phi2)*d2*sbt+wt;

%figure
%plot(t,sat)
%figure
%plot(t,sbt)
%figure
%plot(t,imag(xt))
%figure
%plot(t,abs(xt))

figure
plot(t,real(xt))
xlabel('t')
ylabel('Re(x)')
title('Real part of x')

drf=[0,diff(real(xt))];
figure
plot(t,drf)
line(xlim, [0.5,0.5], 'Color', 'r', 'LineWidth', 0.5); 
line(xlim, [-0.5,-0.5], 'Color', 'r', 'LineWidth', 0.5);
xlabel('t')
ylabel('diff(Re(x))')
title('differentiate of Re(x)')

%%
l=1;
for k=4501:7002
    if abs(drf(k))>0.5
        T(l)=k;
        l=l+1;
    end
end
T

%%
o=1;
for m=2:l-1
    for n=1:l-2
    if T(m)-T(n)==2002
        ta(o)=(T(n)-4501)/ns;
        o=o+1;
    end
    end
end
ta


r=1;
for p=2:l-1
    for q=1:l-2
    if T(p)-T(q)==1001
        tb(r)=(T(q)-4501)/ns;
        r=r+1;
    end
    end
end
tb

taub=min(tb)

for r=1:o-1
if ta(r) ~= taub
    taua=ta(r)
end
end

sa=zeros(1,N);
sa(4501+taua*ns:6502+taua*ns)=1;
a1=4501+taua*ns;
b1=6502+taua*ns;

sb=zeros(1,N);
sb(4501+taub*ns:5501+taub*ns)=1;
sb(5502+taub*ns:6502+taub*ns)=-1;
a2=4501+taub*ns;
b2=6502+taub*ns;

%%
a=min(4501+taua*ns,4501+taub*ns);
b=max(6502+taua*ns,6502+taub*ns);

d1=sum((real(xt(a:b))-(sa(a:b)+sb(a:b))).^2)
d2=sum((real(xt(a:b))-(sa(a:b)-sb(a:b))).^2)
d3=sum((real(xt(a:b))-(-sa(a:b)+sb(a:b))).^2)
d4=sum((real(xt(a:b))-(-sa(a:b)-sb(a:b))).^2)

if min([d1, d2, d3, d4])==d1
    ca=1
    cb=1
    sac=sa;
    sbc=sb;
elseif min([d1, d2, d3, d4])==d2
    ca=1
    cb=0
    sac=sa;
    sbc=-sb;
elseif min([d1, d2, d3, d4])==d3
    ca=0
    cb=1
    sac=-sa;
    sbc=+sb;    
elseif min([d1, d2, d3, d4])==d4
    ca=0
    cb=0
    sac=-sa;
    sbc=-sb;
end

etaa=0.8;
etab=0.8;
phia=-0.8;
phib=-0.8;

ea=0.8;
eb=0.8;
pa=-0.8;
pb=-0.8;

%%
s=sum((abs(xt(a:b)-etaa*exp(i*phia)*sac(a:b)-etaa*exp(i*phia)*sbc(a:b)).^2))

for i1=1:20
    for i2=1:32
        for i3=1:20
            for i4=1:32
                if sum((abs(xt(a:b)-(ea+i1*0.01)*exp(i*(pa+i2*0.05))*sac(a:b)-(eb+i3*0.01)*exp(i*(pb+i4*0.05))*sbc(a:b)).^2))<s
                    s=sum((abs(xt(a:b)-(ea+i1*0.01)*exp(i*(pa+i2*0.05))*sac(a:b)-(eb+i3*0.01)*exp(i*(pb+i4*0.05))*sbc(a:b)).^2));
                    etaa=ea+i1*0.01;
                    phia=pa+i2*0.05;
                    etab=eb+i3*0.01;
                    phib=pb+i4*0.05;
                end
            end
        end
    end
end
   
s
etaa
etab
phia
phib
