Hankel tail

The waveform in equation (31) often arises in practice (as the 2-D Huygens wavelet). Because of the discontinuities on the left side of equation (31), it is not easy to visualize. Thinking again of the time derivative as a convolution with the doublet $(1,-1)/\Delta t$,we imagine the 2-D Huygen's wavelet as a positive impulse followed by negative signal decaying as -t^-3/2. This decaying signal is sometimes called the ``Hankel tail.'' In the frequency domain $-i\omega= \vert\omega \vert e ^ {-i90^\circ}$has a 90 degree phase angle and $\sqrt{-i\omega}= \vert\omega \vert^{1/2} e ^ {-i45^\circ}$has a 45 degree phase angle.

 

# Half order causal derivative.  OK to equiv(xx,yy)
#
subroutine halfdifa( adj, add, n, xx,    yy  )
integer n2, i,       adj, add, n
real omega,                       xx(n), yy(n)
complex cz, cv(4096)
n2=1; while(n2<n) n2=2*n2;         if( n2 > 4096) call erexit('halfdif memory')
do i= 1, n2 {             cv(i) = 0.}
do i= 1, n
        if( adj == 0)   { cv(i) = xx(i)}
        else            { cv(i) = yy(i)}
call adjnull(       adj, add,    xx,n,  yy,n)
call ftu( +1., n2, cv)
                        do i= 1, n2 {
                                omega = (i-1.) * 2.*3.14159265 / n2
                                cz = csqrt( 1. - cexp( cmplx( 0., omega)))
                                if( adj != 0)   cz = conjg( cz)
                                cv(i) = cv(i) * cz
                                }
call ftu( -1., n2, cv)
do i= 1, n
        if( adj == 0)   { yy(i) = yy(i) + cv(i)}
        else            { xx(i) = xx(i) + cv(i)}
return; end

In practice, it is easiest to represent and to apply the 2-D Huygen's wavelet in the frequency domain. Subroutine halfdifa() is provided for that purpose. Instead of using $\sqrt{-i\omega}$ which has a discontinuity at the Nyquist frequency and a noncausal time function, I use the square root of a causal representation of a finite difference, i.e. $\sqrt{1-Z}$,which is well behaved at the Nyquist frequency and has the advantage that the modeling operator is causal (vanishes when t<t₀). Fourier transform is done using subroutine ftu() . Passing an impulse function into subroutine halfdifa() gives the response seen in Figure 7.

 

hankel Figure 7 Impulse response (delayed) of finite difference operator of half order. Twice applying this filter is equivalent to once applying (1,-1).