; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX PRO thm_SCM_casinus_vec, xi,fe,fs,ns,amp,pha,n,nbi, xo, fe_max=fe_max ; ------------------------------------------------------------------ ; * Object : Compute a pure sine from a given signal ; * Class : Data processing for GEOS/UBF ; * Author : P. Robert, CRPE, 1977-1984 ; : Conversion Fortran to IDL, P. Robert, CETP, May 2007 ; : Vectorized by K. Bromund, SPSystems/GSFC, May 2007 ; ------------------------------------------------------------------ ; The input signal contains a high sine signal, with a known ; frequency, superimposed to a useful signal. ; Input : ; xi: array containing the signal ; fe: sampling frequency ; fs: frequency of the sine signal to be removed ; ns: number of spins to use to comput sine wave fit ; Output : ; amp: amplitude of the compute sine signal ; pha: phase (d) of the compute sine signal ; n: number of points in output ; nbi: number of point taken for sine computation ; xo: sine-subtracted signal ; Keyword: ; fe_max: max sampling rate: if fe_max is set, then ; higher sampling rates will be decimated ; to fe_max, and amplitude and phase results will be linearly ; interpolated up to full input rate. ; ------------------------------------------------------------------ ; *** computation of the number of points coresponding to an integer ; number of the sine period to remove n= N_ELEMENTS(xi) if keyword_set(fe_max) && fe gt fe_max then begin fe_priv = fe_max decimation = fe/fe_max if decimation ne fix(decimation) then begin dprint, '*** desinus: warning non-integer decimation' endif n_priv = long(n/decimation) xi_index = lindgen(n_priv)*decimation endif else begin fe_priv = fe n_priv = n xi_index = lindgen(n) endelse nbsp=LONG(FLOAT(n_priv)*fs/fe_priv) ; number of spins in input data nbi =LONG((FLOAT(ns)*fe_priv/fs)+0.5) ;; nuber of points required to ;; fit requested number of spins if not (nbi mod 2) then nbi+=1 ;; algorithm seems to work better w/ odd ; *** test if one have enough points IF (nbsp EQ 0 OR nbi GT n_priv) THEN BEGIN DPRINT, '*** desinus: not enough points to remove the sine' DPRINT, ' at least one sine period is required' DPRINT, ' sampling frequency=',fe_priv DPRINT, ' sine frequency =',fs DPRINT, ' Number of spins requested for fit: ', ns DPRINT, ' number of points required:',nbi, $ ' available=',n_priv amp=0. pha=0. ;Return NaN's for the output values, jmm 23-Oct-2007 xo = xi & xo[*] = !values.f_nan RETURN ENDIF ; compute the sine waves for fitting phase_s = dindgen(nbi)*2*!dpi*fs/fe_priv ;; phase of each sample in kernel ss=sin(phase_s) * hanning(nbi) * 2.0 cc=cos(phase_s) * hanning(nbi) * 2.0 center = nbi/2 ;integer division intentional zs = convol(xi[xi_index], ss) zc = convol(xi[xi_index], cc) amp=2.*sqrt(zs*zs+zc*zc)/FLOAT(nbi) pha=ATAN(zc,zs) ;; this is correct, because we subtract sin rather than ;; cosine wave. pha += phase_s[center] ;; interpolate amplitude and phase to full range and resolution amp = interpol(amp[center:n_priv-center-1], $ xi_index[center:n_priv-center-1], $ lindgen(n)) ;; be careful when interpolating the phase! ;; first, make phase monotonically increasing supplement = 0.d for i = center+1, n_priv-center-1 do begin diffp = pha[i] + supplement - pha[i-1] if diffp gt !dpi then supplement -= 2*!dpi $ else if diffp lt -!dpi then supplement += 2*!dpi pha[i]+=supplement endfor pha = interpol(pha[center:n_priv-center-1], $ xi_index[center:n_priv-center-1], $ lindgen(n)) ; * subtract sine wave while pha is still in pha radians xo = xi - amp * sin(pha) ; * convert phase in degree for output pha *= 180./!dpi pha mod= 360. END ; XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX