;Downsample tplot quantities
;
;
;tname = name of tplot quantity. Can be a size [n] or [n,x]
;sr = new sample rate (S/sec)
;nochange - set if you'd like to overwrite the original tplot variable name
;			Otherwise newname = oldname + '_DS'
;suffix = appended to end of tname. Defaults to 'tname_DS'
;
;win -> apply a hanning window
;ptswin -> number of points on edge of plot for hanning window to go from zero
;			to unity (see hanning_stretch.pro)
;
;Created by Aaron Breneman  10-29-2012
;	Modified: 2012-12-11 : Added anti-aliasing filter
;			  2012-12-14 : Added zero-padding. Greatly improves performance



pro rbsp_downsample,tnames,sr,nochange=nochange,suffix=suffix,win=win,ptswin=ptswin

if ~keyword_set(suffix) then suffix = '_DS'
if ~keyword_set(sr) then sr=1
if ~keyword_set(ptswin) then ptswin = 100.


tn = tnames(tnames)

for j=0,n_elements(tn)-1 do begin


	get_data,tn[j],data=dat,dlimits=dlim,limits=lim
	
	if is_struct(dat) then begin



		;window data
;		if keyword_set(win) then begin
;			han = hanning_stretch(n_elements(dat.x),ptswin)
;;			store_data,'han',data={x:dat.x,y:han}
;;			store_data,'times',data={x:dat.x,y:dat.x-dat.x[0]}
;			dat.y = dat.y * han
;		endif

		;get new times at cadence of "sr"		
		t0 = min(dat.x)
		t1 = max(dat.x)
		newt = dindgen((t1-t0)*sr)/sr + t0
	
		;Test to see if variable is a [n] or [n,x] element array
		tst = size(dat.y,/n_dimensions)
		tst2 = size(dat.y,/dimensions)


		;new Nyquist freq after downsampling will be:
		nyq_new = sr/2.

		;time b/t samples
		dt = dat.x[1]-dat.x[0]


		N = n_elements(dat.x)
		num = lindgen(n_elements(dat.x)/2.)

		;Current range of frequencies in signal is:			
		fcurrent = num/(N*dt)

		;Freqs that will result in aliasing
		badfreqs = where(fcurrent gt nyq_new)
		;Negative freqs that will result in aliasing
		badfreqs_n = reverse(N/2-badfreqs) + N/2


		;Stuff for zero-padding. Greatly speeds up the FFT
		factor = alog(N)/alog(2)
		n_pad_factor = floor(2d^(ceil(factor))) - N
		zero_arr = replicate(0,n_pad_factor)


		if tst eq 1 then begin
		
			;zero pad
			dat_tmp = [dat.y,zero_arr]

		
			;Bandpass original signal to keep freqs only at or below
			;this new Nyquist freq
						
			tmp = fft(dat_tmp)


			;Remove positive freqs that will be aliased
			if badfreqs[0] ne -1 then tmp[badfreqs] = 0.
			
			;Remove negative freqs that will be aliased			
			if badfreqs[0] ne -1 then tmp[badfreqs_n] = 0.


			newsignal = real_part(fft(tmp,1))
			newsignal = newsignal[0:n_elements(newsignal)-n_pad_factor-1]
		
				
			;Interpolate to lower sample rate
			dat2 = interpol(newsignal,dat.x,newt)
		
		
		endif

		


		
		if tst gt 1 then begin					

			dat2 = reform(dblarr(n_elements(newt),tst2[1]))	
			win = hanning(n_elements(dat.x))

			for b=0,tst2[1]-1 do begin


				;zero pad		
				dat_tmp = [dat.y[*,b],zero_arr]


				;Bandpass original signal to keep freqs only at or below
				;this new Nyquist freq
		
				tmp = fft(dat_tmp)
	
	
				;Remove positive freqs that will be aliased
				if badfreqs[0] ne -1 then tmp[badfreqs] = 0.

				;Remove negative freqs that will be aliased			
				if badfreqs[0] ne -1 then tmp[badfreqs_n] = 0.


				newsignal = real_part(fft(tmp,1))
				newsignal = newsignal[0:n_elements(newsignal)-n_pad_factor-1]

				;Interpolate to lower sample rate
				dat2[*,b] = interpol(newsignal,dat.x,newt)
		

			endfor
		endif	
		
		

		
		if ~keyword_set(nochange) then store_data,tn[j]+suffix,data={x:newt,y:dat2},dlimits=dlim,limits=lim $
		else store_data,tn[j],data={x:newt,y:dat2},dlimits=dlim,limits=lim
		
	endif

endfor


end