The helical coordinate

KOLMOGOROFF SPECTRAL FACTORIZATION

Spectral factorization addresses a deep mathematical problem not solved by mathematicians until 1939. Given any spectrum $\vert F(\omega)\vert$ , find a causal time function $f(t)$ with this spectrum. A causal time function is one that vanishes at negative time $t<0$ . We mix spectral factorization with the helix idea to find many applications in geophysical image estimation.

The most abstract method of spectral factorization is of the Russian mathematician A.N.Kolmogoroff. I include it here, because it is by far the fastest, so much so that giant problems become practical, such as the solar physics example coming up.

Given that $C(\omega)$ Fourier transforms to a causal function of time, it is next proven that $e^C$ Fourier transforms to a causal function of time. Its filter inverse is $e^{-C}$ . Grab yourself a cup of coffee and hide yourself away in a quiet place while you focus on the proof in the next paragraph.

A causal function $c_\tau$ vanishes at negative $\tau$ . Its $Z$ transform $C(Z) = c_0 + c_1 Z + c_2 Z^2 + c_3 Z^3 +\cdots$ , with $Z=e^{i\omega\Delta t}$ is really a Fourier sum. Its square $C(Z)^2$ convolves a causal with itself so it is causal. Each power of $C(Z)$ is causal, therefore, $e^C=1+C+C^2/2+\cdots$ , a sum of causals, is causal. The time-domain coefficients for $e^C$ could be computed putting polynomials into power series or computed faster with Fourier transforms (by understanding $C(Z=e^{i\omega\Delta t})$ as an FT.) By the same reasoning, the wavelet $e^C$ has inverse $e^{-C}$ which is also causal. A causal with a causal inverse is said to be ``minimum phase.'' The filter $1-Z/2$ with inverse $1+Z/2+Z^2/4+\cdots$ is minimum phase because both are causal, and they multiply to make the impulse ``1'', so are mutually inverse. The delay filter $Z^5$ has the noncausal inverse $Z^{-5}$ which is not causal (output before input).

The next paragraph defines ``Kolmogoroff spectral factorization.'' It arises in applications where one begins with an energy spectrum $\vert r\vert^2$ and factors it into an $r e^{i\phi}$ times its conjugate. The inverse Fourier transform of that $r e^{i\phi}$ is causal.

Relate amplitude $r=r(\omega)$ and phase $\phi=\phi(\omega)$ to a causal time function $c_\tau$ .

$$\vert r\vert e^{i\phi} = e^{\ln\vert r\vert}e^{i\phi} \ =\ e^{\ln\,\vert r\vert + i\phi} \ =\ e^{c_0+c_1Z+c_2Z^2+c_3Z^3+\cdots} \ =\ e^{\sum_{\tau=0} c_\tau Z^\tau}$$ (7)

Given a spectrum $r(\omega)$ , we find a filter with that spectrum. Because $r(\omega)$ is a real even function of $\omega$ , so is its logarithm. Let the inverse Fourier transform of $\ln \vert r(\omega)\vert$ be $u_\tau$ , where $u_\tau$ is a real even function of time. Imagine a real odd function of time $v_\tau$ .

$$\vert r\vert e^{i\phi} = e^{\ln\,\vert r\vert + i\phi} \ =\ e^{\sum_\tau (u_\tau+v_\tau) Z^\tau}$$ (8)

The phase $\phi(\omega)$ transforms to $v_\tau$ . We can assert causality by choosing $v_\tau$ so that $u_\tau+v_\tau=0$ $(=c_\tau)$ for all negative $\tau$ . This choice defines $v_\tau$ at negative $\tau$ . Since $v_\tau$ is odd, it is also known at positive lags. More simply, $v_\tau$ is created when $u_\tau$ is multiplied by a step function of size 2. This causal exponent $(c_0,c_1,\cdots)$ creates a causal filter $\vert r\vert e^{i\phi}$ with the specified spectrum $r(\omega)$ .

We easily manufacture an inverse filter by changing the polarity of the $c_\tau$ . This filter is also causal by the same reasoning. Thus, these filters are causal with a causal inverse. Such filters are commonly called ``minimum phase.''

Spectral factorization arises in a variety of contexts. Here is one: Rain drops showering on a tin roof create for you a signal with a spectrum you can compute, but what would be the sound of a single drop, the wavelet of a single drop? Spectral factorization gives the answer. Divide this wavelet out from the data to get a record of impulses, one for each rain drop (theoretically!). Similarly, the boiling surface of the sun is coming soon.

Everyone has their own favorite Fourier transform code, so why am I offering mine? Because you MUST get the scale factors correct. Few worries if you accidentally replace $e^C$ by $2e^C$ , because your humble plotting program might do that. But, if you accidentally replace $e^C$ by $e^{2C}$ , you have squared it!

subroutine ftu( signi, nx, cx )	          # Fourier transform
#   complex Fourier transform with traditional scaling (FGDP)
#
#               1         nx          signi*2*pi*i*(j-1)*(k-1)/nx
#   cx(k)  =  -------- * sum cx(j) * e
#              scale     j=1             for k=1,2,...,nx=2**integer
#
#  scale=1 for forward transform signi=1, otherwise scale=1/nx
integer nx, i, j, k, m, istep
real    signi, arg
complex cx(nx), cmplx, cw, cdel, ct
i=1;  while( i<nx) i=2*i
if( i != nx )    call erexit('ftu: nx not a power of 2')
do i= 1, nx
        if( signi<0.)
                cx(i) = cx(i) / nx
j = 1;  k = 1
do i= 1, nx {
        if (i<=j) { ct = cx(j); cx(j) = cx(i); cx(i) = ct }
        m = nx/2
        while (j>m && m>1) { j = j-m; m = m/2 }         # "&&" means .AND.
        j = j+m
        }
repeat {
        istep = 2*k;   cw = 1.;   arg = signi*3.14159265/k
        cdel = cmplx( cos(arg), sin(arg))
        do m= 1, k {
                do i= m, nx, istep
                        { ct=cw*cx(i+k);  cx(i+k)=cx(i)-ct;  cx(i)=cx(i)+ct }
                cw = cw * cdel
                }
        k = istep
        if(k>=nx) break
        }
return; end

The ftu fast Fourier transform code has a restriction that the data length must be a power of 2. Zero time and frequency are the first point in the vector, then positive times, and then negative times.

It is an exercise for the student to show that a complex-valued time function has a positive spectrum that is nonsymmetrical in frequency, but it may be factored with the same code.

The helical coordinate