Doxygen Source Code Documentation
thd_detrend.c File Reference
#include "mrilib.h"Go to the source code of this file.
Functions | |
| void | THD_const_detrend (int npt, float *xx, float *xx0) |
| void | get_linear_trend (int npt, float *xx, float *f0, float *f1) |
| void | THD_linear_detrend (int npt, float *far, float *xx0, float *xx1) |
| void | get_quadratic_trend (int npt, float *xx, float *f0, float *f1, float *f2) |
| void | THD_quadratic_detrend (int npt, float *far, float *xx0, float *xx1, float *xx2) |
| void | THD_cubic_detrend (int npt, float *far) |
| void | THD_normalize (int npt, float *far) |
| void | THD_generic_detrend (int npt, float *far, int polort, int nort, float **ort) |
Function Documentation
|
||||||||||||||||||||
|
Definition at line 27 of file thd_detrend.c. References x0. Referenced by absfft_func(), linear_filter_trend(), and THD_linear_detrend().
00028 {
00029 double t1,t3,t10 , x0,x1 ;
00030 int ii ;
00031
00032 if( npt < 2 || xx == NULL || f0 == NULL || f1 == NULL ) return ;
00033
00034 x0 = xx[0] ; x1 = 0.0 ;
00035 for( ii=1 ; ii < npt ; ii++ ){
00036 x0 += xx[ii] ;
00037 x1 += xx[ii] * ii ;
00038 }
00039
00040 t1 = npt*x0; t3 = 1.0/npt; t10 = npt*npt;
00041
00042 *f0 = (float)(2.0/(npt+1.0)*t3*(2.0*t1-3.0*x1-x0));
00043 *f1 = (float)(-6.0/(t10-1.0)*t3*(-x0-2.0*x1+t1));
00044 return ;
00045 }
|
|
||||||||||||||||||||||||
|
Definition at line 76 of file thd_detrend.c. Referenced by THD_quadratic_detrend().
00077 {
00078 double x0,x1,x2 , N=npt ;
00079 int ii ;
00080
00081 if( npt < 3 || xx == NULL || f0 == NULL || f1 == NULL || f2 == NULL ) return;
00082
00083 x0 = xx[0] ; x1 = x2 = 0.0 ;
00084 for( ii=1 ; ii < npt ; ii++ ){
00085 x0 += xx[ii] ;
00086 x1 += (xx[ii] * ii) ;
00087 x2 += (xx[ii] * ii) * ii ;
00088 }
00089
00090 *f0 = ( 3.0*(3.0*N*N-3.0*N+2.0) * x0
00091 -18.0*(-1.0+2.0*N) * x1
00092 +30.0 * x2 ) / (N*(N+2.0)*(N+1.0)) ;
00093
00094 *f1 = ( -18.0*(-1.0+2.0*N) * x0
00095 +12.0*(-1.0+2.0*N)*(8.0*N-11.0) * x1 /((N-1.0)*(N-2.0))
00096 -180.0 * x2 /(N-2.0) )
00097 / (N*(N+2.0)*(N+1.0)) ;
00098
00099 *f2 = ( 30.0 * x0
00100 -180.0 * x1 / (N-2.0)
00101 +180.0 * x2 / ((N-1.0)*(N-2.0)) ) / (N*(N+2.0)*(N+1.0)) ;
00102 return ;
00103 }
|
|
||||||||||||||||
|
Definition at line 11 of file thd_detrend.c.
00012 {
00013 int ii ; float xbar ;
00014
00015 if( npt < 2 || xx == NULL ) return ;
00016
00017 xbar = 0.0 ;
00018 for( ii=0 ; ii < npt ; ii++ ) xbar += xx[ii] ;
00019 xbar /= npt ;
00020 for( ii=0 ; ii < npt ; ii++ ) xx[ii] -= xbar ;
00021 if( xx0 != NULL ) *xx0 = xbar ;
00022 return ;
00023 }
|
|
||||||||||||
|
Cubic detrend a float array in place. Definition at line 134 of file thd_detrend.c. References far.
00135 {
00136 register int ii ;
00137 float g0,g1,g2,g3 , f0,f1,f2,f3 , t1,t2,t5,t8 , t95,t56,t22,t25,txx ;
00138
00139 if( npt < 5 || far == NULL ) return ;
00140
00141 t8 = npt*npt ; t2 = npt-1.0 ; t5 = t2*(npt-2.0) ;
00142 t95 = 0.05*t5*(npt-3.0) ;
00143 t56 = 0.16666667*t5 ;
00144 t22 = 0.5*t2 ;
00145 t25 = 1.5*t2 ;
00146 txx = 0.6*t8-1.5*npt+1.1 ;
00147
00148 g0=g1=g2=g3=0.0 ;
00149 for( ii=0 ; ii < npt ; ii++ ){
00150 t1 = ii*ii ;
00151 f1 = ii - t22 ;
00152 f2 = t1 - t2*ii + t56 ;
00153 f3 = t1*(ii - t25) + txx*ii - t95 ;
00154
00155 g0 += far[ii] ;
00156 g1 += far[ii] * f1 ;
00157 g2 += far[ii] * f2 ;
00158 g3 += far[ii] * f3 ;
00159 }
00160 g0 *= (1.0/npt) ;
00161 g1 *= 12.0/(npt*(t8-1.0)) ;
00162 g2 *= 180.0/(npt*(t8-1.0)*(t8-4.0)) ;
00163 g3 *= 2800.0/(npt*(t8-1.0)*(t8-4.0)*(t8-9.0)) ;
00164
00165 for( ii=0 ; ii < npt ; ii++ ){
00166 t1 = ii*ii ;
00167 f1 = ii- t22 ;
00168 f2 = t1 - t2*ii + t56 ;
00169 f3 = t1*(ii - t25) + txx*ii - t95 ;
00170
00171 far[ii] -= ( g0 + g1*f1 + g2*f2 + g3*f3 ) ;
00172 }
00173
00174 return ;
00175 }
|
|
||||||||||||||||||||||||
|
Detrend a vector with a given polort level, plus some others. ------------------------------------------------------------------------- Definition at line 200 of file thd_detrend.c. References far, fit, free, lsqfit(), malloc, and ref. Referenced by main().
00202 {
00203 int ii,jj , nref ;
00204 float **ref , *fit , xmid , xfac , val ;
00205
00206 /* check inputs */
00207
00208 if( npt <= 1 || far == NULL ) return ;
00209 if( nort > 0 ){
00210 if( ort == NULL ) return ;
00211 for( jj=0 ; jj < nort ; jj++ ) if( ort[jj] == NULL ) return ;
00212 }
00213 if( polort < 0 ) polort = -1 ;
00214 if( nort < 0 ) nort = 0 ;
00215
00216 nref = polort+1+nort ;
00217 if( nref == 0 || nref >= npt-1 ) return ;
00218
00219 ref = (float **) malloc( sizeof(float *) * nref ) ;
00220 xmid = 0.5*(npt-1) ; xfac = 1.0 / xmid ;
00221 for( jj=0 ; jj <= polort ; jj++ ){
00222 ref[jj] = (float *) malloc( sizeof(float) * npt ) ;
00223 switch( jj ){
00224 case 0:
00225 for( ii=0 ; ii < npt ; ii++ ) ref[jj][ii] = 1.0 ;
00226 break ;
00227
00228 case 1:
00229 for( ii=0 ; ii < npt ; ii++ ) ref[jj][ii] = xfac*(ii-xmid) ;
00230 break ;
00231
00232 case 2:
00233 for( ii=0 ; ii < npt ; ii++ ){
00234 val = xfac*(ii-xmid) ; ref[jj][ii] = val*val ;
00235 }
00236 break ;
00237
00238 case 3:
00239 for( ii=0 ; ii < npt ; ii++ ){
00240 val = xfac*(ii-xmid) ; ref[jj][ii] = val*val*val ;
00241 }
00242 break ;
00243
00244 default:
00245 for( ii=0 ; ii < npt ; ii++ ){
00246 val = xfac*(ii-xmid) ; ref[jj][ii] = pow(val,(double)(jj)) ;
00247 }
00248 break ;
00249 }
00250 }
00251 for( jj=0 ; jj < nort ; jj++ ) /* user supplied refs */
00252 ref[polort+1+jj] = ort[jj] ;
00253
00254 fit = lsqfit( npt , far , NULL , nref , ref ) ;
00255
00256 if( fit != NULL ){ /* good */
00257 for( ii=0 ; ii < npt ; ii++ ){
00258 val = far[ii] ;
00259 for( jj=0 ; jj < nref ; jj++ )
00260 val -= fit[jj] * ref[jj][ii] ;
00261 far[ii] = val ;
00262 }
00263 free(fit) ;
00264 } else {
00265 fprintf(stderr,"THD_generic_detrend: lsqfit fails - no detrending!\n") ;
00266 }
00267
00268 for( jj=0 ; jj <= polort ; jj++ ) free(ref[jj]) ;
00269 free(ref) ; return ;
00270 }
|
|
||||||||||||||||||||
|
Definition at line 53 of file thd_detrend.c. References far, and get_linear_trend(). Referenced by main(), normalize_floatvector(), and THD_dataset_tshift().
00054 {
00055 register int ii ;
00056 float f0 , f1 ;
00057
00058 if( npt < 3 || far == NULL ) return ;
00059
00060 get_linear_trend( npt , far , &f0 , &f1 ) ;
00061
00062 far[0] -= f0 ;
00063 for( ii=1 ; ii < npt ; ii++ ) far[ii] -= (f0 + f1*ii) ;
00064
00065 if( xx0 != NULL ) *xx0 = f0 ;
00066 if( xx1 != NULL ) *xx1 = f1 ;
00067
00068 return ;
00069 }
|
|
||||||||||||
|
Definition at line 181 of file thd_detrend.c. References far. Referenced by main().
00182 {
00183 register int ii ;
00184 register float fac ;
00185
00186 if( npt <= 0 || far == NULL ) return ;
00187
00188 fac = 0 ;
00189 for( ii=0 ; ii < npt ; ii++ ) fac += far[ii]*far[ii] ;
00190 if( fac == 0.0 ) return ;
00191 fac = 1.0 / sqrt(fac) ;
00192 for( ii=0 ; ii < npt ; ii++ ) far[ii] /= fac ;
00193 return ;
00194 }
|
|
||||||||||||||||||||||||
|
Definition at line 111 of file thd_detrend.c. References far, and get_quadratic_trend(). Referenced by main().
00113 {
00114 register int ii ;
00115 float f0 , f1 , f2 ;
00116
00117 if( npt < 4 || far == NULL ) return ;
00118
00119 get_quadratic_trend( npt , far , &f0 , &f1 , &f2 ) ;
00120
00121 far[0] -= f0 ;
00122 for( ii=1 ; ii < npt ; ii++ ) far[ii] -= ( (f2*ii + f1)*ii + f0 ) ;
00123
00124 if( xx0 != NULL ) *xx0 = f0 ;
00125 if( xx1 != NULL ) *xx1 = f1 ;
00126 if( xx2 != NULL ) *xx2 = f2 ;
00127
00128 return ;
00129 }
|
