PRO Conefilt, mov, dx, dt, vel, HIGHPASS=highpass, ADD_APO=adap ;+ ; NAME: ; CONEFILT ; PURPOSE: ; k-w filtering of a 2D time series (in memory) ; CALLING SEQUENCE: ; Conefilt, data, dx, dt, Vel ; INPUTS: ; Data: 3-D datacube (x,y,t) with images. Will be ; overwritten by result! ; dx, dt: (spatial and temporal) stepwidth of data, eg 51, 30 ; vel: cutoff Velocity in km/s ; OPTIONAL PARAMETERS: ; ; KEYWORDS: ; ; OUTPUTS: ; ; RESTRICTIONS: ; ; PROCEDURE: ; ; MODIFICATION HISTORY: ; 20-Oct-1999 P.Suetterlin, SIU ; 26-Oct-1999 First and last image are somehow "smeared". Try ; to put an artificial interpolated image before ; and after the series. ; 08-Feb-2001 Rewrite: more general, fix bugs ; 16-Apr-2004 make addition of interpolated images optional: ; they influence the powerspectra of the filtered ; cube... ;- on_error, 2 IF n_params() LT 4 THEN BEGIN message, "Use: Conefilt, data, dx, dt, Vel", /informal return ENDIF s = size(mov) IF s(0) NE 3 THEN BEGIN message, "Input data must be a 3d data cube (x,y,t)!" return ENDIF sx = s(1) sy = s(2) num = s(3) dtype = s(4) memuse = (sx*sy*(num+2))*8 ;;; stop if the data cube needs more than so many MB ;RR changed 1500 to 2000 = 2/3 of my 3GB Tosh memory IF memuse GT 2000*1024l^2 THEN BEGIN message, "Data cube exceeds 2 GB - use conefilt_file instead!" return ENDIF ;;; create an array with the radial wavenumbers. ;;; horizontal and vertical binsize are assumed to be equal. ;;; code taken from IDL routine dist. kr = fltarr(sx, sy, /nozero) k_ny = 0.5/dx x = findgen(sx) x = (x < (sx-x)) x = ((x/max(x))*k_ny)^2 sy2 = float(sy/2) tmp = k_ny/sy2 FOR i=0., sy2 DO BEGIN y = sqrt(x + (tmp*i)^2) kr(0, i) = y IF i NE 0 THEN kr(0, sy-i) = y ENDFOR ;;; temporal boundary conditions: Insert an additional (identical) ;;; frame at the beginning and end ;;; ?? Or use apodisation ?? ;;; No, checked apodisation - does not work as it leaves artifacts ;;; in the data. The added frames however give problems if you ;;; look at the powerspectra of the cube later: The masked-out bins ;;; are no longer zero, most likely due to numerical problems with ;;; the discrete frequency bins that are different. ;;; This 'apodisation' is optional now. IF keyword_set(adap) THEN BEGIN sm1 = (2.*mov(*, *, 0) + mov(*, *, num-1))/3. sm2 = (mov(*, *, 0) + 2.*mov(*, *, num-1))/3. mov = [[[sm1]], [[temporary(mov)]], [[sm2]]] n = num+2 ENDIF ELSE n = num ;;; Do the FFT mov = fft(temporary(mov), -1) ;;; temporal frequency points n2 = n/2 f = findgen(n2+1)/n2*(0.5/dt) ;;; take out the cone with kr < f/Vel FOR i=1, n/2 DO BEGIN IF keyword_set(highpass) THEN $ mask = kr LE (f(i)/Vel) $ ELSE $ mask = kr GT (f(i)/Vel) mov(*, *, i) = mov(*, *, i) * mask mov(*, *, n-i) = mov(*, *, n-i) * mask ENDFOR ;;; FFT back, convert to float and cut off first and last image mov = fft(mov, 1, /over) mov = float(temporary(mov)) IF keyword_set(adap) THEN mov = (temporary(mov))(*, *, 1:n-2) ;;; Convert to original data type if not float IF dtype EQ 1 THEN mov = byte(mov < 255 > 0) IF dtype EQ 2 THEN mov = fix(mov < 32767 > (-32767)) END