; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
; deconvolib
;
; Library of routines for deconvolution of UBF waveforms
;
; Allows despinning and calibration in nT of waveforms to low
; frequency (typically 0.1-10 Hz)
;
; Assumes existence of function which returns the
; complex gain of the antennas in V/nT for a given frequency,
; antenna,calibration file for a given satellite, and sample
; frequency. If given a negative frequency, it must return the
; complex conjugate of the gain. The form of the function must be
; gainant(f, ix, isat, fe),
;
; Based on deconvolib by P. Robert CNRS/CETP
;
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
pro thm_scm_apotrap, s, ns
; thm_scm_apotrap: apply a trapeziodal window to a waveform. waveform
; may be real or complex
; Based on P. Robert's apotrap in deconvolib.f90
na = ns/16
ap = dindgen(na)/double(na-1)
ap_ind = lindgen(na)
s[ap_ind] *= ap
ap_ind = ns -1 - ap_ind
s[ap_ind] *= ap
end
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
pro thm_scm_corgain, s, ns, df, fe, gainant, ix, pathfil,frq, init=init
; thm_scm_corgain: apply correction for antenna gain to complex
; spectrum s
; Based on P. Robert's corgain in deconvolib.f90
if keyword_set(init) then s = dcomplexarr(ns)+complex(1.0,0)
frq = findgen(ns)*df
frq[ns/2+1:*] -= float(ns)*df
;; IDL version of gainant handles negative frequencies.
if keyword_set(gainant) then $
c = call_function(gainant, frq, ix, pathfil, fe) $
else c = complex(1.0, 0)
smallc = where(abs(c) lt 1.e-6, nsmallc)
if nsmallc gt 0 then c[smallc] = complex(1.e-6, 0.)
s /= c
;; FFT of real function must be real at f=0
s[0]=abs(s[0])*((real_part(s[0]) lt 0) ? -1 : 1)
end
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
pro thm_scm_coretar, s, ns, fe, tau
; thm_scm_coretar: apply correction for sample delay of tau seconds
; to complex spectrum s
; Based on P. Robert's coretar in deconvolib.f90
t = double(ns)/fe
ns2=ns/2
n = lindgen(ns)
n[ns2+1:*] -= ns
teta = 2.0d*!dpi*n*double(tau)/t
srot = exp(dcomplex(0, 1)*teta)
s *= srot
end
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
pro thm_scm_filtspe, s, ns, df, f1, f2
; thm_scm_filtspe: filter a complex spectrum s below f1 and above f2
; Based on P. Robert's filtspe in deconvolib.f90
f = findgen(ns)*df
f[ns/2+1:*] -= float(ns)*df
f = abs(f)
; light smoothing below f1 [sic]
;arg = (f/(2*df))**2
;arg = arg < 37.
;af2 = 1.-exp(-arg)
; changed January 21, 2002: rectangular filter
af2 = 1.
filt = where(f lt f1 or f gt f2, nfilt)
if nfilt gt 0 then s[filt] = complex(0., 0.)
;filt = where(f ge f1 && f le f2, nfilt)
;if nfilt gt 0 then s[filt] *= af2
end
; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX
pro thm_scm_deconvo_vec, xio, nbp, nk, f1, f2, fe, gainant, ix, pathfil, tau, $
edge_truncate=edget, edge_wrap=edgew, edge_zero=edgez, $
plotk=plotk
;+
;Name:
; thm_scm_deconvo_vec
;Purpose:
; continuous calibration accomplished by
; convolveing a signal with a kernel.
; The kernel is derived by taking the fft of the
; inverse of the frequency reponse. The kernel is also constructed
; to account for any sample delay.
;Inputs:
; xio: input/output waveform
; nbp: number of points in xio
; nk: number of points in kernel: power of 2.
; f1: low frequency
; f2: high frequency
; fe: sample frequency
;gainant: gain function, such as thm_scm_gainant
; ix: antenna number (1,2,3)
;pathfil: calibration file name
; tau: time correction
;-
cx = ['','X','Y','Z']
df = fe/float(nk)
;; pre-centering
xavg = total(xio)/nbp
xio -= xavg
;; generate kernel that compensates for gain and sample delay.
;; first get desired complex frequency reponse of filter
s = dcomplexarr(nk)+dcomplex(1.0, 0)
;; corgain: inverse of antenna gain (s,nk,df,fe,gainant,ix,pathfil)
thm_scm_corgain, s, nk, df, fe, gainant, ix, pathfil, frq
;; coretar: correct for sample delay (s,nk,fe,tau)
thm_scm_coretar, s, nk, fe, tau
;; filtspe: bandpass filter between f1 and f2 (s,nbp,df,f1,f2)
thm_scm_filtspe, s, nk, df, f1, f2
;; inverse fft to obtain kernel
ks = fft(s, 1)
kernel = real_part(ks)
;; if we did everything right, the imaginary part truly is negligible
pr = total(kernel*kernel)
pi = total(imaginary(ks)*imaginary(ks))
if pi/pr gt 1.e-5 then begin
print, '*** thm_scm_deconvo: Imag/Real for impulse reponse for ', cx[ix], $
' =', pi/pr
endif
;; zero phase of the kernel is at index 0. shift that to center, nk/2,
;; for use with IDL convol routine.
kernel = shift(kernel, nk/2)
;; As this is a continuous calibration, the window must be applied to the
;; kernel, rather than to the waveform.
kernel_orig = kernel
kernel *= hanning(nk, /double)
;; note: application of window introduces a slight offset, which must be
;; removed from the signal afterwards.
;; Correcting for the offset in the kernel itself
;; would nullify the benefit the of window.
if keyword_set(plotk) then begin
;; visualize effects of window on frequency response:
;; (we expect to see an effect on the phase, when compared
;; to ideal correction spectrum, due to centering of kernel)
!p.multi=[0,0,4]
;; plot kernel - compare windowed to unwindowed
plot, kernel_orig, /xstyle, color=1
oplot, kernel, color=2
re = max(kernel)
oplot, [nk/2, nk/2], [-re, +re], color=3
;; plot amplitude response, compare shifted windowed and unwindowed to ideal.
nk21 = nk/2+1
plot, shift(frq,-nk21), shift(abs(s),-nk21), /ylog, $
yrange=[0.01,max(abs(s))], /xstyle
oplot, shift(frq,-nk21), shift(abs(fft(kernel_orig)),-nk21), color=1
oplot, shift(frq,-nk21), shift(abs(fft(kernel)),-nk21), color=2
;; plot amplitude response, compare shifted windowed and unwindowed to ideal.
plot, shift(frq,-nk21), shift(atan(s,/ph)*!radeg,-nk21), /xstyle
oplot, shift(frq,-nk21),shift(atan(fft(kernel_orig),/ph)*!radeg,-nk21),color=1
oplot, shift(frq,-nk21), shift(atan(fft(kernel),/ph)*!radeg,-nk21), color=2
;; plot group delay
plot, shift(frq,-nk21), -(shift(atan(s, /ph),-nk21-1) - $
shift(atan(s,/ph),-nk21))/(df*2*!dpi), /xstyle
endif
;; perform the convolution.
;; notes on edge behavior:
;; default: zero output when kernel overlaps edge
;;; /edge_zero: usually good -- can emphasize low frequency trends, i.e.
;;; artifiacts of despin
;;; /edge_wrap; similar to edge zero for cal signal.
;;; /edge_truncate: usually bad
kernel /= nk
xio = convol(xio, kernel, /center, $
edge_t=edget, edge_w=edgew, edge_z=edgez)
;; post-centering of waveform (necessary because of window)
xavg = total(xio)/nbp
xio -= xavg
end