module cross_wilson{ # Wilson's factorization of cross-correlation use helicon use polydiv integer, parameter, private :: n = 1024 real, dimension( 2*n-1), private :: cross, ab, cc real, dimension( n), private :: b, c contains subroutine cross_factor( niter, s0, ss, a0, aa, b0, bb) { integer, intent( in) :: niter # Newton iterations real, intent( in) :: s0 # in: cross-correlat. type( filter), intent( in) :: ss # real, intent( in out) :: a0, b0 # out: min. phase type( filter), intent( out) :: aa, bb # real :: eps integer :: i, stat cross = 0.; cross( n) = s0; b( 1) = 1. # initialize cross( n+ss%lag) = ss%flt # cross-correlation do i = 1, niter { call polydiv_init( 2*n-1, aa) # divide polynoms stat = polydiv_lop( .false., .false., cross, ab) # ab = S/A # ab = ab/(a0*b0) # scale call polydiv_init( 2*n-1, bb) # divide polynoms stat = polydiv_lop( .true., .false., cc, ab) # cc = S/(AB') b( 2:n) = cc( n+1:2*n-1)/cc( n) # b = + side # b( 1) = 0.5*(cc( n) + 1.) # of (1 + cc) eps = maxval(abs(b(2:n))) call helicon_init( aa) # mutliply polynoms stat = helicon_lop( .false., .false., b, c) # c = A b # c = c*a0; a0 = c( 1); aa%flt = c(1+aa%lag)/a0 # scale aa%flt = c(1+aa%lag) b( 2:n) = cc( n-1:1:-1)/cc( n) # b = - side eps = max(maxval(abs(b(2:n))),eps) write( 0,*) i, eps # "L1 norm" if( eps < epsilon(a0)) exit # convergence call helicon_init( bb) # mutliply polynoms stat = helicon_lop( .false., .false., b, c) # c = B b # c = c*b0; b0 = c( 1); bb%flt = c(1+bb%lag)/b0 # scale bb%flt = c(1+bb%lag) } call polydiv_close() } }