Source Code for the S Transform of a Real Function

; Returns the complex localization matrix for positive frequencies

; The users passes the time series, and the width factor of the localizing window

; ie the window = k periods (k=1 suggested)

; result = local(timeseries,k)

• FUNCTION local, time_series , factor
• width=1.0 ; declares a floating point variable
• temp = size(time_series)
• length = temp(1) ; finds the length of the time series array
• b = complexarr(length) ; allocates a complex array with 'length' elements
• gw = complexarr(length) ; allocates a complex array with 'length' elements
• spe_length = length/2
• loc = complexarr(length,spe_length+1) ; allocates a complex two dimensional array ; with 'length' by 'spe_length+1' elements
• h = fft(time_series,-1) ; multi-radix fft routine defined below
• loc(*,0) = 0 ; zero frequency is set to zero (subtract average)
• for f = 1.0, spe_length do begin
• width = factor * length/f ; width of window in the time domain
• gw = g_window(length,width) ; user defined function described below
• h = shift(h,-1) ; shifts elements in an array
• b = h * gw ; multiplies the spectrum by the 'voice gaussian'
• loc(*,f) = fft(b,1) ; a voice of the localization matrix (complex)
• endfor
• return,loc
• end

; Gaussian Window Function called by local()

; Accepts variables length and width

; Returns the spectrum of the Gaussian window

• FUNCTION g_window, len, w
• if w ne 0.0 then sigma = len/(2*!PI*w) $ ;sigma = standard deviation of voice gaussian
• else print,'w is zero!'
• g = fltarr(len) ; allocate memory for voice gaussian
• compg = complex(len) ;allocates a complex array of length 'len'
• for i = 0.0,(len-1) do begin
• exponent = (i-(len/2))^2/(2*sigma^2) ; calculate exponent of gaussian
• if exponent lt 25 then g(i)=exp(-exponent) $ ; if exponent large then e^(-exp)=0
• else g(i)=0.0 ; { usually not necessary }
• endfor
• g = shift(g,-len/2) ; positive half first, then negative half
• return,complex(g,0) ; returns the complex-valued voice gaussian
• end

; Inverse of the local() function

; Returns the original time series (as a complex array)

; Accepts the positive frequency half of the complex S transform result

; ie this accepts the output of the local() function

; Usage: inverse_result = ilocal(s-matrix)

; for example, in the following, inverse_result will be the original time series

; IDL> result = local(timeseries,1)

; IDL> inverse_result = ilocal(result)

; then REAL[inverse_result] == timeseries

• FUNCTION ilocal, loc
• tempsize = size(loc) ; get size of the complex array 'loc'
• rlength = tempsize(1) ; length of row
• clength = tempsize(2) ; length of column
• lsum = complexarr(clength)
• for i = 0,clength-1 do lsum(i) = total( (temp=loc(*,i)) ) ; averages the voices to give the ;Fourier spectrum (+ freq only!)
• tempsize=0 & temp=0
• lspe = complexarr(rlength)
• for i = 0,clength-1 do lspe(i) = lsum(i)
• for i = clength,rlength-1 do
• lspe(i) =complex(float(lsum(2*clength-i-2)),-imaginary(lsum(2*clength-i-2)) )
• ; mirrors the negative frequencies (real even,imag odd)
• ; thus the array 'lspe' is the fourier transform of ;the array "timeseries"
• return, fft(lspe/rlength,1) ;inverse fft to return the time series (complex valued)
• ; note, the imaginary part of the time series will be ;identically equal to zero.
• end

IDL Routines

fft(array,+-1)

{for spectrum = fft( array,-1) }

and the inverse Fourier transform is:

{ array = fft(spectrum,1) }

Both functions return a complex array of length N. It is stored in the usual way of positive frequencies in 0-N/2 followed by the negative frequencies.

SHIFT(array, n)