;+
;FUNCTION: THM_LSP_REMOVE_SPIKES, t, xx, per, nwin=nwin, nfit=nfit
;
; NOT FOR GENERAL USE. CALLED BY THM_EFI_REMOVE_SPIKES
; ONLY FOR ISOLATED PARTICLE OR WAVE BURSTS
; NOT FOR ENTIRE ORBIT.
;
;PURPOSE:
; Remove the non-physical spiky signals in the efw data.
;
;INPUT:
; t -NEEDED. Time array
; x -NEEDED. Data
; per -NEEDED. Spin period.
;
;KEYWORD:
; nwin -OPTIONAL. Number of points in spike search window. DFLT = 16
; spikesig -OPTIONAL. Sigma of spikes. DFLT = 5
; sigmin -OPTIONAL. Minimun of sigma. DFLT = 0 mV/m.
; nfit -OPTIONAL. Number of points in the fit window. DFLT = 16
;
;HISTORY:
; 2008-11-08: Created by Jianbao Tao at LASP@CU-Boulder.
; 2008-11-10: The head comment was added.
; 2009-03-30: REE. Rewrite to test epoch analysis.
; 2009-04-30: REE. Rewrite to use median smooth.
;-
function thm_lsp_remove_spikes, t, xx, per, nwin=nwin, spikesig=spikesig, $
sigmin=sigmin, nfit=nfit, talk=talk, diagnose=diagnose
; CHECK KEYWORDS
if not keyword_set(nwin) then nwin = 16L
if not keyword_set(spikesig) then spikesig = 5.0d
if n_elements(sigmin) EQ 0 then sigmin = 0.01d
if not keyword_set(nfit) then nfit = 16L
x = xx
; CALCULATE AVERAGE SIGNAL OVER PERIOD
tt = t-t(0)
dt = t(1:*)-t(*)
mdt = median(dt)
phs = tt/per
npps = long(round(2.0d*per/mdt)) ; SET NUMBER OF POINTS PER SPIN AT ~2X DT
nspins = ceil(max(phs)) ; NUMBER OF SPINS
ntot = long(nspins) * long(npps) ; NUMBER OF PHASE POINTS
phx = dindgen(ntot)/npps
xp = interpol(x, phs, phx)
xr = reform(xp, npps, nspins)
xe = total(xr,2)/ nspins
if keyword_set(talk) then $
plot, phx(0:npps-1), xe, xtit = 'Spin Phase (/2pi)', $
title = 'Spin-Averaged Value', ytit = 'Amplitude'
; FIND SPIKES
xel = [xe(npps-nwin/2:npps-1), xe, xe(0:nwin/2-1)]
xel2 = abs(xel-thm_lsp_median_smooth(xel,nwin))
sig = abs(xel2)/(stddev(xel2(nwin/2: npps-nwin/2-1))+sigmin)
sdum = sig ; ISOLATE MAXIMUMS IN SIG
FOR i= nwin/2-1, npps-nwin/2 DO BEGIN
smax = max(sdum(i-nwin/2+1:i+nwin/2))
if sdum(i) EQ smax then sig(i)=smax else sig(i) = 0
ENDFOR
sig = sig(nwin/2: npps-nwin/2-1)
ppk = where(sig GE spikesig, nppk)
; CHECK PPK FOR ERRORS
IF nppk EQ 0 then BEGIN
print, 'NO SPIKES FOUND'
IF keyword_set(talk) then BEGIN
xyouts, 0.25, 0.75, 'NO SPIKES FOUND', charsize=2, /normal
wait, 1.0
ENDIF
return, x
ENDIF
IF nppk GT 4 THEN BEGIN
print, 'TOO MANY SPIKES - CANNOT REMOVE'
return, x
ENDIF
; PEAKS
tpks = ppk*per/double(npps)
ipks = round(tpks/dt)
iper = float(per/mdt)
IF keyword_set(talk) then BEGIN
oplot, phx(ppk), xe(ppk), psym=6, col = 250
xyouts, 0.25, 0.75, 'SPIKES', charsize=2, /normal
wait, 0.25
ENDIF
; CALCULATE TIMES / INDECIES
tel = [phx(npps-nwin/2:npps-1)-1.d, phx(0: npps+nwin/2-1)]
nfit = 16
good = 0
bad = 0
; REMOVE PEAKS
FOR j = 0, n_elements(ipks)-1 DO BEGIN
; DO A PRELIMINARY FIT TO DETERMINE LIMITS
xwin = xel[ppk(j)-nfit/2+nwin/2:ppk(j)+nfit/2+nwin/2]
twin = tel[ppk(j)-nfit/2+nwin/2:ppk(j)+nfit/2+nwin/2]
fit = gaussfit(twin, xwin, A, nterms=6, chisq=chisq, sigma=sigma)
; DIAGNOSTIC PLOT
IF keyword_Set(diagnose) then BEGIN
plot, twin, xwin
oplot, twin, A(0)*exp( -((twin-A(1))/A(2))^2 / 2.0d) + $
poly(twin, A(3:5)) , col = 2
ENDIF
A0 = A(0)
A1 = A(1)
A0min = 0
A0max = 0
if A0 GT 0 then A0max = 5.0*A0 else A0min = 5.0*A0
istart = float((ipks(j)-nfit/2) > 0)
iend = round(istart) + nfit
ispin = 0
; NEXT LOOP
WHILE iend LT n_elements(x)-1 DO BEGIN
xwin = x(round(istart):iend)
twin = tt(round(istart):iend) - tt(round(istart))
fit = gaussfit(twin, xwin, A, nterms=6, chisq=chisq, sigma=sigma)
cent = round(A(1)/mdt) - nfit/2
; DIAGNOSTIC PLOT
if keyword_set(diagnose) then plot, twin, xwin
; FIND GOOD FITS AND REMOVE ELSE SUBTRACT TYPICAL
IF (chisq LE 3) AND (abs(cent) LE 2) AND $
(A(0) GT A0min) AND (A(0) LT A0max) then BEGIN
xwin = xwin - A(0)*exp( -((twin-A(1))/A(2))^2 / 2.0d)
good = good + 1
; DIAGNOSTIC PLOT
IF keyword_Set(diagnose) then BEGIN
oplot, twin, A(0)*exp( -((twin-A(1))/A(2))^2 / 2.0d) + poly(twin, A(3:5)) , col = 2
xyouts, 0.25, 0.75, 'GOOD FIT', charsize=2, /normal
ENDIF
ENDIF else BEGIN ; NOTCH
A = ladfit(twin, xwin)
xwin(nwin/4:3*nwin/4) = poly(twin(nwin/4:3*nwin/4), A)
bad = bad + 1
; DIAGNOSTIC PLOT
IF keyword_Set(diagnose) then BEGIN
oplot, twin, poly(twin, A) , col = 2
xyouts, 0.25, 0.75, 'NOTCH', charsize=2, /normal
ENDIF
ENDELSE
; DIAGNOSTIC PLOT
IF keyword_Set(diagnose) then BEGIN
oplot, twin, xwin, col=6
wait, 0.25
ENDIF
; RESET X
x(round(istart):iend) = xwin
istart = istart+iper; + (( (round(A(1)/mdt) - nfit/2) > (-1)) < 1); REMOVED
iend = istart + nfit
ispin = ispin + 1
ENDWHILE
END
print, good, ' SPIKES REMOVED (GOOD FITS)', bad, ' SPIKES NOTCHED (FAILED FITS)'
return, x
end
; DIAGNOSTICS
;plot, twin, xwin
;oplot, twin, A(0)*exp( -((twin-A(1))/A(2))^2 / 2.0d) + $
;A(3) + A(4)*twin + A(5)*twin*twin, col=2
;AAV1 = A0(1)*per + ispin*per - tt(round(istart))
;oplot, twin, AAV(0)*exp( -((twin-AAV1)/AAV(2))^2 / 2.0d) + $
;A(3) + A(4)*twin + A(5)*twin*twin, col=4
;oplot, twin, xwin, col=6
;print, AAV1, A(1)
;if abs(cent) GT 2 then print, 'POOR CENTER'
;if chisq GT 3 then print, 'POOR FIT'
;if A(0) LT A0min or A(0) GT A0max then print, 'POOR AMPLITUDE'
; OLD NOTCHER
;great = 0
;AAV = fltarr(6)
; FIND GREAT FITS TO DETERMINE AAV
;IF (chisq LT 1) AND (abs(cent) LT 2) AND $
;(A(0) GT A0min) AND (A(0) LT A0max) then BEGIN
;AAV = (AAV*great + A)/(1+great)
;great=great+1
;ENDIF
;ENDIF else BEGIN
;AAV1 = A0(1)*per + ispin*per - tt(round(istart))
;xwin = xwin - AAV(0)*exp( -((twin-AAV1)/AAV(2))^2 / 2.0d)
;bad = bad + 1
;ENDELSE
; OLD PEAK FINDER
; FIND SPIKES
;slide = nwin/2
;xel = [xe,xe]
;ppk = 0
;FOR i = 0L, npps-1, slide DO BEGIN
; wstart = i
; wend = wstart + nwin
; win = xel[wstart: wend]
; xmed = median(win)
; xstd = stddev(win)
; xpk = max(abs(win-xmed), ipeak)
; print, xpk, xstd, xpk/xstd
; if (xpk GT (xstd*3)) AND (ipeak GT nwin/4) AND (ipeak LE 3*nwin/4) then $
; ppk = [ppk, wstart+ipeak]
;ENDFOR
; CHECK PPK FOR ERRORS
;IF n_elements(ppk) EQ 1 then BEGIN
; print, 'NO SPIKES FOUND'
; return, x
;ENDIF
;IF n_elements(ppk) GT 5 THEN BEGIN
; print, 'TOO MANY SPIKES - CANNOT REMOVE'
; return, x
;ENDIF
; PEAKS
;ppk = ppk(1:*)