;+
; NAME:
; THM_LSP_GET_SPEC (PROCEDURE)
;
; PURPOSE:
; Get the power spectrum density (PSD) of a given tplot variable.
;
; CALLING SEQUENCE:
; thm_lsp_get_spec, tvar, units = units, prefix = prefix, fftlen = fftlen, $
; yrange = yrange
;
; ARGUMENTS:
; tvar: (INPUT, REQUIRED) The name of a tplot variable to calculate the
; PSD from.
;
; KEYWORDS:
; units: (INPUT, OPTIONAL) A string of the units of the data in the tvar. By
; default, it is obtained from dlim.data_att.units, or 'unknown' if
; dlim.data_att.units is not available.
; prefix: (INPUT, OPTIONAL) The prefix for tplot variables that contain the
; resulting PSD. By default, prefix = tvar.
; fftlen: (INPUT, OPTIONAL) The number of data points to FFT. By default,
; fftlen = 512.
; yrange: (INPUT, OPTIONAL) The yrange for the spectra. By default, yrange =
; [1, Nyquist] in units of Hz.
;
; EXAMPLES:
; See thm_crib_lsp_get_spec.pro which is typically located under
; TDAS_DIR/idl/themis/examples.
;
; HISTORY:
; 2010-03-07: Created by JBT, CU/LASP.
;
;-
pro thm_lsp_get_spec, tvar, units = units, prefix = prefix, fftlen = fftlen, $
yrange = yrange
; Check tvar.
thm_lsp_checkarg_tvar, tvar, error, errmsg
if error ne 0 then begin
print, 'THM_LSP_GET_PSD: ' + errmsg
print, 'Exiting...'
return
endif
get_data, tvar, data=data, dlim = dlim
; Check units.
if not keyword_set(units) then begin
str_element, dlim, 'data_att', success = s0
if not s0 then units = 'unknown' else begin
str_element, dlim.data_att, 'units', success = s1
if not s1 then units = 'unknown' else units = dlim.data_att.units
endelse
endif
tmp1 = size(units, /type) ne 7
tmp2 = n_elements(units) ne 1
tmp = tmp1 + tmp2
if tmp ne 0 then begin
print, 'THM_LSP_GET_PSD: ' + $
'Valid units of the input tplot variable must be a string scalar. '+$
'Exiting...'
return
endif
units0 = units[0]
; Check prefix.
if n_elements(prefix) eq 0 then pre = tvar else begin
pre = prefix
if size(pre, /type) ne 7 then pre = tvar
if n_elements(pre) ne 1 then pre = tvar
if pre ne prefix then begin
print, 'THM_LSP_GET_PSD: ' + $
'WARNING -- No valid prefix is given. The default prefix is used.'
endif
endelse
pre = pre[0]
; Check fftlen.
wrongfftlen = 0
if n_elements(fftlen) ne 1 then begin
wrongfftlen = 1
fftlen = 512
endif
if size(fftlen, /type) lt 2 or size(fftlen, /type) gt 5 then begin
wrongfftlen = 1
fftlen = 512
endif
if wrongfftlen then $
print, 'THM_LSP_GET_PSD: ' + $
'WARNING -- No valid fftlen is given. The default fftlen (512) is used.'
dt = median(data.x[1:*] - data.x)
tarr = data.x
dim = size(data.y,/dim)
nyquist = 0.5 / dt
if not keyword_set(yrange) then yrange = [1, nyquist]
; Generate spectrum names
ncomp = dim[1]
if ncomp eq 3 then name = pre + ['xspec', 'yspec', 'zspec'] else begin
tmp = indgen(ncomp) + 1
name = pre + string(tmp, for='(I0)') + 'spec'
endelse
; Some constants.
nslide = fftlen/2 ; Number of points to slide to next section.
time_units = 1. ; [sec]
freq_units = 1. ; [Hz]
data_units = 1. ; whatever units of the input is.
; Determine the frequencies of the power spectrum.
df = 1. / (fftlen * dt)
nf = fftlen/2 + 1
farr = dindgen(nf) * df
fftdata_units = data_units / sqrt(freq_units) ; units of fft(data).
psd_units = fftdata_units^2 ; PSD units.
; Generate a Hann window for fft.
win = hanning(fftlen, /double) ; The Hann window.
btrange = thm_jbt_get_btrange(tvar, nb=nb, tind=tind)
; Get the time stamps and frequency bins.
for ib = 0L, nb -1L do begin
ista = tind[ib, 0]
iend = tind[ib, 1]
nt = iend - ista + 1L
if nt lt fftlen then continue
if n_elements(ibstart) eq 0 then ibstart = ib
nsec = long((nt - fftlen) / nslide) + 1L
tmptime = dindgen(nsec) * (dt * nslide) + tarr[nslide+ista]
tmpfreq = transpose(rebin(farr, nf, nsec))
if n_elements(time) lt 1 then time = tmptime else time=[time,tmptime]
if n_elements(freq) lt 1 then freq = tmpfreq else freq=[freq, tmpfreq]
endfor
if n_elements(time) lt 1 then begin
print, 'THM_LSP_GET_SPEC: ' + $
'No continuous section is longer than a fftlen. No spectrum is stored.'
print, 'Exiting...'
return
endif
psdunits = units0 + '!E2!N / Hz'
coord = cotrans_get_coord(tvar)
att = {units:psdunits, coord_sys:coord, data_type:'calibrated'}
; dlim = {spec:1B, log:1B, data_att:att, ylog:1, zlog:1, $
; ytitle:xpsdname, ztitle:psdunits, yrange:yrange}
for icom = 0L, ncomp-1L do begin
psdname = name[icom]
dlim = {spec:1B, log:1B, data_att:att, ylog:1, zlog:1, $
ztitle:psdunits, yrange:yrange}
for ib = 0L, nb-1L do begin
ista = tind[ib, 0]
iend = tind[ib, 1]
nt = iend - ista + 1L
if nt lt fftlen then continue
nsec = long((nt - fftlen) / nslide) + 1L
dat = data.y[ista:iend, icom]
tmppsd = dblarr(nf, nsec)
for isec = 0L, nsec - 1 do begin
; print, ''
; print, '# ', isec + 1, ' out of ', nsec
iista = isec * nslide
iiend = iista + fftlen - 1L
xdat = dat[iista:iiend]
fx = fft(xdat * win)
;wx- energy conservation factor for xdat
wx = total(xdat^2) / total((xdat*win)^2)
sqrfx = real_part(fx * conj(fx))
; Get PSD of the two series.
tmppsd[0,isec] = sqrfx[0]
tmppsd[1:fftlen/2-1,isec] = sqrfx[1:fftlen/2-1] + $
reverse(sqrfx[fftlen/2+1:fftlen-1])
tmppsd[fftlen/2,isec] = sqrfx[fftlen/2]
tmppsd[*,isec] = tmppsd[*,isec] * wx
endfor
; Average tmppsd.
xpsd = tmppsd
if nsec eq 1 then xpsd = transpose(xpsd) else begin
if nsec eq 2 then begin
xpsd[*,0] = (tmppsd[*,0] + tmppsd[*,1]*0.5) / 1.5
xpsd[*,nsec-1] = (tmppsd[*,nsec-1] + tmppsd[*,nsec-2]*0.5) / 1.5
xpsd = transpose(xpsd)
endif else begin
xpsd[*,0] = (tmppsd[*,0] + tmppsd[*,1]*0.5) / 1.5
xpsd[*,nsec-1] = (tmppsd[*,nsec-1] + tmppsd[*,nsec-2]*0.5) / 1.5
for isec = 1L, nsec - 2 do begin
xpsd[*,isec] = (total(tmppsd[*,isec-1:isec+1], 2) + $
tmppsd[*,isec])/4.
endfor
xpsd = transpose(xpsd)
endelse
endelse
if ib eq ibstart then psd=xpsd else psd=[psd, xpsd]
endfor
psd = psd * psd_units
store_data, psdname, data={x:time, y:psd, v:freq}, dlim=dlim
endfor
end