;+ ; 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