module kolmogroff_mod
use scaletrace_mod
use myfft_mod
use ftu_mod
implicit none

include "fftw3.f"
integer*8       ,save		:: planF,planB
real   , allocatable            :: in(:)
complex, allocatable            :: out(:)
complex, allocatable            :: in2(:)
integer		,save		:: n
real		,save		:: scaleF, scaleB

contains
!+=============
subroutine mpwave_init(cx)
	complex			:: cx(:)
	n =size(cx)
	scaleF=sqrt( 1./n)
	scaleB=sqrt( 1./n)
	allocate(in(n),in2(n),out(n/2+1))
	call sfftw_plan_dft_r2c_1d(planF,n,in,out,FFTW_ESTIMATE)
	call sfftw_plan_dft_c2r_1d(planB,n,out,in,FFTW_ESTIMATE)
end subroutine
!+=============
subroutine mpwave(cx)       		!# minimum phase equivalent wavelet
	integer i            		!# input:  cx = any wavelet
	complex cx(:)                   !# output: cx = min phase wavelet 
	in = cx				! FFT weighting
	call sfftw_execute(planF); 	out=out*scaleF
	out=out*conjg(out)
	call kolmogoroff()
	cx=in
	return; 
end subroutine
!+=============
subroutine mpwaveA(cx)       		!# spectrum
	integer i            		!# input:  cx = any wavelet
	complex cx(:)                   !# output: cx = min phase wavelet 
	in = cx		 		! FFT weighting
	call sfftw_execute(planF); out=out*scaleF
	call kolmogoroff()
	cx=in
	return; 
return; 
end subroutine
!+=============
subroutine kolmogoroff()  		!# Spectral factorization.
	integer i       		!# input:  cx = spectrum !# output: cx = min phase wavelet
	where(out/=0.)out = clog(out); 
	call sfftw_execute(planB);	in =in *scaleB
	in(1    ) = in(1    ) / 2.
	in(1+n/2) = in(1+n/2) / 2.
	in(1+n/2+1:n) = 0.
	call sfftw_execute(planF); 	out = out*scaleF
	out = cexp(out) 
	call sfftw_execute(planB); 	 in = in*scaleB 
	in =in/in(1)*1.791
	return; 
end subroutine

subroutine mpwave_close()
	deallocate(in,in2,out)
	call sfftw_destroy_plan(planF);	
	call sfftw_destroy_plan(planB);
end subroutine
end module



