c *****************************************
c  DCMP-SD algorithm (2 elements, 2 waves)
c *****************************************
      parameter (ndata = 500)
      complex*16 xin(2,0:ndata),y(0:ndata),ww(2,0:ndata+1),wtmp(2)
      complex*16 dwave, uwave, dsig, usig, zj
      complex*16 C(2),P(2,2),F(2), H
      real*8 pow1, pow2, ang1, ang2, ang_s, pi, derad
      real*8 sinr(0:ndata), cal_sinr, mu, noise
      real 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 = 1.0d0
      noise = 0.01d0
      st = sqrt(noise/2.0)
      ang1 = 30.0d0
      ang2 = -60.0d0

      open( 9,file= 'sinr-dcmp-sd0002.dat')
      open(10,file='weight-dcmp-sd0002.dat')
      write(10,*) ndata+1

c constraint vector & constraint response
      ang_s = 30.0d0
      C(1) = 1.0d0
      C(2) = cdexp(-zj*pi*sin(ang_s*derad))
      H = (1.0d0,0.0d0)

c Projection matrix and F vector
      P(1,1) = 1.0d0 - 0.5d0*C(1)*dconjg(C(1))
      P(1,2) = - 0.5d0*C(1)*dconjg(C(2))
      P(2,2) = 1.0d0 - 0.5d0*C(2)*dconjg(C(2))
      P(2,1) = dconjg(P(1,2))
      F(1) = 0.5d0 * dconjg(H) * C(1)
      F(2) = 0.5d0 * dconjg(H) * C(2)

c  initial state & stepsize
      ww(1,0) = F(1)
      ww(2,0) = F(2)
      mu = 0.003d0


c  weight update loop
      do i = 0, ndata
         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) = 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) = dsig*cdexp(-zj*pi*sin(ang1*derad))
     &        + usig*cdexp(-zj*pi*sin(ang2*derad))
     &        + dble(thmr) + zj*dble(thmi)

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

c  SINR
         sinr(i) = cal_sinr(ww(1,i))

c  weight update
         wtmp(1) = ww(1,i) - mu*xin(1,i)*dconjg(y(i))
         wtmp(2) = ww(2,i) - mu*xin(2,i)*dconjg(y(i))
         ww(1,i+1) = P(1,1)*wtmp(1)+P(1,2)*wtmp(2) + F(1)
         ww(2,i+1) = P(2,1)*wtmp(1)+P(2,2)*wtmp(2) + F(2)

      end do
      
c  results
      do i = 0, ndata
         write( 9,*) i, float(sinr(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  output SINR
      real*8 function cal_sinr(weight)
      complex*16 weight(*), y_s, y_u
      complex*16 V_s(2), V_u(2), zj/(0.0d0,1.0d0)/
      real*8 pow1, pow2, ang1, ang2, pi, derad, noise, Pnout
      common /rdata/pow1,ang1,pow2,ang2,noise
      pi = acos(-1.0d0)
      derad = pi/180.0d0
      V_s(1) = 1.0d0
      V_s(2) = cdexp(-zj*pi*sin(ang1*derad))
      V_u(1) = 1.0d0
      V_u(2) = cdexp(-zj*pi*sin(ang2*derad))
      y_s = V_s(1)*dconjg(weight(1))+V_s(2)*dconjg(weight(2))
      y_u = V_u(1)*dconjg(weight(1))+V_u(2)*dconjg(weight(2))
      Pnout = noise *
     &     (weight(1)*dconjg(weight(1))+weight(2)*dconjg(weight(2)))
      cal_sinr = pow1*abs(y_s)**2 / (pow2*abs(y_u)**2 + Pnout)
      cal_sinr = 10.0d0*log10(cal_sinr)
      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