PRO vikfft,pathnum,orbit,id,f,bsfft,wpeaks,_EXTRA=e

bs=1
rdimg2,pathnum,orbit,id,h,i0,i1,stat,stream=bs

IF stat EQ 1 THEN BEGIN
    isize=SIZE(i1)
    nimg=2
    ENDIF ELSE BEGIN
    rdimg,pathnum,orbit,id,h,i
    isize=SIZE(i)
    IF isize(N_ELEMENTS(isize)-1) EQ 0 THEN BEGIN
	MESSAGE,'Cannot find the desired image in the specified path',/INFORMATIONAL
	RETURN
	ENDIF ELSE BEGIN
	bs=REFORM(i,isize(1)*isize(2))
        nimg=1
	ENDELSE
    ENDELSE

pixperrow=isize(1)

n=N_ELEMENTS(bs)
t=1./32768.
n21=n/2+1
f=LINDGEN(n)
f(n21)=n21-n+FINDGEN(n21-2)
f=f/(n*t)

IF !D.NAME EQ 'X' THEN BEGIN
    wwidth=MAX([2*nimg*pixperrow,640])
    WINDOW,2,XSIZE=640,YSIZE=720
    ENDIF ELSE DEVICE,YOFFSET=1,YSIZE=9,/INCHES
!P.MULTI=[0,0,3,0,0]

PLOT,bs
multi2=!P.MULTI

bsfft=FFT(bs,-1)
absfft=ABS(bsfft)
mabsfft=MEDIAN(absfft,n/256)

IF !D.NAME EQ 'X' THEN BEGIN
    wwidth=MAX([2*nimg*pixperrow,640])
    WINDOW,0,XSIZE=wwidth,YSIZE=960
    ENDIF ELSE DEVICE,YOFFSET=1,YSIZE=9,/INCHES
;!P.POSITION=[0.15,0.05,0.95,0.65]
!P.MULTI=[0,0,4,0,0]

tstring=STRING(orbit,id,FORMAT='(I0,1X,F6.3)')
!X.TITLE='Frequency (Hz)'
!P.TITLE=tstring
PLOT_IO,SHIFT(f,-n21),SHIFT(absfft,-n21),_EXTRA=e $
       ,YTITLE='Power spectral density (units?)'
PLOT_IO,SHIFT(f,-n21),SHIFT(mabsfft,-n21),_EXTRA=e $
       ,YTITLE='Median-smoothed PSD'

wpeaks=WHERE(absfft/mabsfft GT 5)
rowpers=32768./(nimg*pixperrow)
cycperrow=f(wpeaks)/rowpers
PLOT,HISTOGRAM(ABS(cycperrow),MIN=0,BIN=1),PSYM=10 $
    ,XTITLE='Cycles per image row',XRANGE=[0,100]

PRINT,'Enter frequencies (in cycles/row) to be removed'
rms=''
READ,rms
nrmf=nwrds(rms)
rmf=FLTARR(nrmf)
READS,rms,rmf

fbsfft=bsfft
mr=MEDIAN(FLOAT(bsfft))
mi=MEDIAN(IMAGINARY(bsfft))

FOR j=0,nrmf-1 DO BEGIN
    wrmp=WHERE(ABS(ABS(f/rowpers)-rmf(j)) LT 0.5)
    fbsfft(wrmp)=REPLICATE(COMPLEX(mr,mi),N_ELEMENTS(wrmp))
    ENDFOR

PLOT_IO,SHIFT(f,-n21),SHIFT(ABS(fbsfft),-n21),_EXTRA=e $
       ,YTITLE='Filtered PSD',NOERASE=0
fbs=FFT(fbsfft,1)

multi0=!P.MULTI
WSET,2
!P.MULTI=multi2
PLOT,FLOAT(fbs)
PLOT,IMAGINARY(fbs)
multi2=!P.MULTI
WSET,0
!P.MULTI=multi0

STOP

!ORDER=1
IF nimg EQ 1 THEN mtv,i,0 ELSE BEGIN
    mtv,i0,0
    mtv,i1,1
    ENDELSE

IF nimg EQ 1 THEN BEGIN
    fi=REFORM(fbs,isize(1),isize(2))
    fi=BYTE(ROUND(ABS(fi) < 255))
    mtv,fi,1
    ENDIF ELSE BEGIN
    fi=REFORM(fbs,2,isize(1)*isize(2))
    fi1=REFORM(fi(0,*),isize(1),isize(2))
    fi0=REFORM(fi(1,*),isize(1),isize(2))
    fi1=BYTE(ROUND(ABS(fi1) < 255))
    fi0=BYTE(ROUND(ABS(fi0) < 255))
    mtv,fi0,2
    mtv,fi1,3
    ENDELSE

RETURN
END