c*********************************
c     1-2 LS-CMA algorithm
c*********************************
      parameter (ndata = 100, indat = 20)
      complex*16 rxx(2,2),rxxi(2,2)
      complex*16 xin(2,ndata,indat),y(ndata,indat),ww(2,0:ndata),
     &     dr(ndata,indat)
      complex*16 dwave, uwave, dsig, usig, zj, er
      complex*16 rtmp(2),update(2)
      real*8 pow1, pow2, ang1, ang2, pi, derad
      real*8 ee(0:ndata), alpha
      real noise, st, thmr, thmi
      integer i, idst/584287/, iust/350107/, inst/455782/
      common /rdata/pow1,ang1,pow2,ang2,noise
      pi = acos(-1.0d0)
      derad = pi / 180.0d0
      zj = (0.0d0, 1.0d0)

c  signal parameters
      pow1 = 1.0d0
      pow2 = 0.5d0
      noise = 0.01d0
      st = sqrt(noise/2.0)
      ang1 = 30.0d0
      ang2 = -60.0d0

      open( 9,file='error-lscma.dat')
      open(10,file='weight-lscma.dat')
      write(10,*) ndata

c initial state
      ww(1,0) = (1.0d0, 0.0d0)
      ww(2,0) = (0.0d0, 0.0d0)
c
      alpha = 0.0d0
      sigma = 1.0d0
c
c  weight update loop
      do i = 1, ndata
         rxx(1,1) = (0.0d0,0.0d0)
         rxx(2,2) = (0.0d0,0.0d0)
         rxx(1,2) = (0.0d0,0.0d0)
         rxx(2,1) = (0.0d0,0.0d0)
         rtmp(1) = (0.0d0,0.0d0)
         rtmp(2) = (0.0d0,0.0d0)

c  loop for update vector calculation
         do in = 1, indat
            dsig = sqrt(pow1)*dwave(idst)
            usig = sqrt(pow2)*uwave(iust)

c  input at 1st element
            call normal(inst,0.0,st,thmr)
            call normal(inst,0.0,st,thmi)
            xin(1,i,in) = dsig + usig + dble(thmr)+zj*dble(thmi)

c  input at 2nd element
            call normal(inst,0.0,st,thmr)
            call normal(inst,0.0,st,thmi)
            xin(2,i,in) = dsig*cdexp(-zj*pi*sin(ang1*derad))
     &           + usig*cdexp(-zj*pi*sin(ang2*derad))
     &           + dble(thmr) + zj*dble(thmi)

c  array output
            y(i,in) = dconjg(ww(1,i-1))*xin(1,i,in)
     &           +dconjg(ww(2,i-1))*xin(2,i,in)


c  reference signal
            if ( abs(y(i,in)) .lt. 1.0d-5 ) then
               dr(i,in) = sigma
            else
               dr(i,in) = sigma / abs(y(i,in)) * y(i,in)
            end if

c  correlation matrix
            rxx(1,1)=rxx(1,1)
     &           + xin(1,i,in)*dconjg(xin(1,i,in))
            rxx(1,2)=rxx(1,2)
     &           + xin(1,i,in)*dconjg(xin(2,i,in))
            rxx(2,1)=rxx(2,1)
     &           + xin(2,i,in)*dconjg(xin(1,i,in))
            rxx(2,2)=rxx(2,2)
     &           + xin(2,i,in)*dconjg(xin(2,i,in))

c  error signal
            er = dr(i,in) - y(i,in) 

            rtmp(1) = rtmp(1) + xin(1,i,in)*dconjg(er)
            rtmp(2) = rtmp(2) + xin(2,i,in)*dconjg(er)

         end do

c  averaged cost function
         ee(i-1) = 0.0d0
         do in = 1, indat
            ee(i-1) = ee(i-1)+(abs(y(i,indat))**2 - sigma**2)**2
         end do
         ee(i-1) = ee(i-1) / dble(indat)

c  calculation of Rxx^{-1}
         rxxd = (rxx(1,1)+alpha)*(rxx(2,2)+alpha) - rxx(1,2)*rxx(2,1)
         rxxi(1,1) = (rxx(2,2)+alpha) / rxxd
         rxxi(2,2) = (rxx(1,1)+alpha) / rxxd
         rxxi(1,2) = -rxx(1,2) / rxxd
         rxxi(2,1) = -rxx(2,1) / rxxd

         update(1)=rxxi(1,1)*rtmp(1)+rxxi(1,2)*rtmp(2)
         update(2)=rxxi(2,1)*rtmp(1)+rxxi(2,2)*rtmp(2)

c  weight update
         ww(1,i) = ww(1,i-1) + update(1)
         ww(2,i) = ww(2,i-1) + update(2)

      end do

c  results
      do i = 0, ndata-1
         write( 9,*) i, float(ee(i))
         write(10,*) ww(1,i),ww(2,i)
      end do

      close(9)
      close(10)

      stop
      end

c  desied signal
      complex*16 function dwave(idst)
      complex*16 zj/(0.0d0,1.0d0)/
      real*8 pi,phase
      real pp
      integer idst
      pi = acos(-1.0d0)
      call random(idst,pp)
      phase = 2.0d0*pi*pp
      dwave = cdexp(zj*phase)
      return
      end

c  undesired signal
      complex*16 function uwave(iust)
      complex*16 zj/(0.0d0,1.0d0)/
      real*8 pi,phase
      real pp
      integer iust
      pi = acos(-1.0d0)
      call random(iust,pp)
      phase = 2.0d0*pi*pp
      uwave = cdexp(zj*phase)
      return
      end

c  normal random
      subroutine normal(x,av,st,rn)
      integer x
      real av,st,rn,r1,r2,w1,w2
      call random(x,r1)
      call random(x,r2)
      w1=log(r1)
      w1=sqrt(-2.*w1)
      w2=sin(6.2831852*r2)
      rn=av+st*w1*w2
      return
      end

c  uniform random
      subroutine random(x,ru)
      integer x
      real ru
      x=x*65539
      if (x.lt.0) x=(x+2147483647)+1
      ru=x*0.4656613e-9
      return
      end