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
!------------
integer, parameter   				:: FFTW_R2HC=0,FTW_HC2R=1,FFTW_DHT=2,FFTW_REDFT00=3,FFTW_REDFT01=4
integer, parameter   				:: FFTW_REDFT10=5,FFTW_REDFT11=6,FFTW_RODFT00=7,FFTW_RODFT01=8,FFTW_RODFT10=9
integer, parameter   				:: FFTW_RODFT11=10,FFTW_FORWARD=-1,FFTW_BACKWARD=+1,FFTW_DESTROY_INPUT=1
integer, parameter   				:: FFTW_UNALIGNED=2,FFTW_CONSERVE_MEMORY=4,FFTW_EXHAUSTIVE=8,FTW_PRESERVE_INPUT=16
integer, parameter   				:: FFTW_PATIENT=32,FFTW_ESTIMATE=64,FFTW_ESTIMATE_PATIENT=128
integer, parameter   				:: FFTW_BELIEVE_PCOST=256,FFTW_DFT_R2HC_ICKY=512,FFTW_NONTHREADED_ICKY=1024
integer, parameter   				:: FFTW_NO_BUFFERING=2048,FFTW_NO_INDIRECT_OP=4096,FFTW_ALLOW_LARGE_GENERIC=8192
integer, parameter   				:: FFTW_NO_RANK_SPLITS=16384,FFTW_NO_VRANK_SPLITS=32768,FFTW_NO_VRECURSE=65536
!------------
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



