00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 #include "mrilib.h"
00008 
00009 /******************* Routines to shift two rows at a time ********************/
00010 
00011 /*---------------------------------------------------------------------------
00012    Set the interpolation method for shifting:
00013    input is one of MRI_NN, MRI_LINEAR, MRI_CUBIC, or MRI_FOURIER.
00014 -----------------------------------------------------------------------------*/
00015 
00016 typedef void (*shift_func)(int,int,float,float *,float,float *) ;
00017 static  shift_func shifter      = fft_shift2 ;
00018 static  int        shift_method = MRI_FOURIER ;
00019 
00020 void SHIFT_set_method( int mode )
00021 {
00022    shift_method = mode ;
00023    switch( mode ){
00024       default:          shift_method = MRI_FOURIER ;  /* fall thru */
00025       case MRI_FOURIER: shifter = fft_shift2   ; break ;
00026 
00027       case MRI_LINEAR:  shifter = lin_shift2   ; break ;
00028       case MRI_CUBIC:   shifter = cub_shift2   ; break ;
00029       case MRI_QUINTIC: shifter = quint_shift2 ; break ;  /* Nov 1998 */
00030       case MRI_HEPTIC:  shifter = hept_shift2  ; break ;  /* Nov 1998 */
00031 
00032       case MRI_NN:      shifter = nn_shift2    ; break ;  /* experimental */
00033       case MRI_TSSHIFT: shifter = ts_shift2    ; break ;  /* Dec 1999 */
00034    }
00035    return ;
00036 }
00037 
00038 int SHIFT_get_method( void ){ return shift_method ; }
00039 
00040 /*--------------------------------------------------------------------------
00041    The main entry point - note that g can be NULL, but f cannot.
00042 ---------------------------------------------------------------------------*/
00043 
00044 void SHIFT_two_rows( int n, int nup, float af, float * f, float ag, float * g )
00045 {
00046    shifter( n,nup,af,f,ag,g ) ; return ;
00047 }
00048 
00049 /*--------------------------------------------------------------------------
00050    Shift 2 rows at a time with the FFT:
00051      n   = length of a row
00052      nup = length to use for FFTs (must be even)
00053      af  = shift for row f
00054      ag  = shift for row g
00055    Input and output arrays are f[n] and g[n].  (Note: g may be NULL.)
00056 ----------------------------------------------------------------------------*/
00057 
00058 #define ZFILL
00059 #define RECUR
00060 
00061 void fft_shift2( int n, int nup, float af, float * f, float ag, float * g )
00062 {
00063    static int nupold=0 , nuptop=0 ;
00064    static complex * row=NULL , * cf=NULL , * cg=NULL ;
00065 
00066    int ii , nby2=nup/2 , n21=nby2+1 ;
00067    complex fac , gac ;
00068    float sf , sg , dk ;
00069 #ifdef RECUR
00070    complex csf , csg ;
00071 #endif
00072 
00073 ENTRY("fft_shift2") ;
00074 
00075    /* 15 Mar 2001: shift too big ==> return all zeros */
00076 
00077    if( (af < -n || af > n) && (ag < -n || ag > n) ){
00078       for( ii=0 ; ii < n ; ii++ ) f[ii] = g[ii] = 0.0 ;
00079       EXRETURN ;
00080    }
00081 
00082    /* make new memory for row storage? */
00083 
00084    if( nup > nuptop ){
00085       if( row != NULL ){ free(row) ; free(cf) ; free(cg) ; }
00086       row = (complex *) malloc( sizeof(complex) * nup ) ;
00087       cf  = (complex *) malloc( sizeof(complex) * n21 ) ;
00088       cg  = (complex *) malloc( sizeof(complex) * n21 ) ;
00089       nuptop = nup ;
00090    }
00091 
00092    /* FFT the pair of rows */
00093 
00094    if( g != NULL )
00095       for( ii=0 ; ii < n ; ii++ ){ row[ii].r = f[ii] ; row[ii].i = g[ii] ; }
00096    else
00097       for( ii=0 ; ii < n ; ii++ ){ row[ii].r = f[ii] ; row[ii].i = 0 ; }
00098 
00099 #ifdef ZFILL
00100    for( ii=n ; ii < nup ; ii++ ){ row[ii].r = row[ii].i = 0.0 ; }
00101 #else
00102    if( nup > n ){
00103       sf = 0.5 * (row[0].r + row[n-1].r) ; sg = 0.5 * (row[0].i + row[n-1].i) ;
00104       for( ii=n ; ii < nup ; ii++ ){ row[ii].r = sf ; row[ii].i = sg ; }
00105    }
00106 #endif
00107 
00108    csfft_cox( -1 , nup , row ) ;
00109 
00110    /* untangle FFT coefficients from row into cf,cg */
00111 
00112    cf[0].r = 2.0 * row[0].r ; cf[0].i = 0.0 ;  /* twice too big */
00113    cg[0].r = 2.0 * row[0].i ; cg[0].i = 0.0 ;
00114    for( ii=1 ; ii < nby2 ; ii++ ){
00115       cf[ii].r =  row[ii].r + row[nup-ii].r ;
00116       cf[ii].i =  row[ii].i - row[nup-ii].i ;
00117       cg[ii].r =  row[ii].i + row[nup-ii].i ;
00118       cg[ii].i = -row[ii].r + row[nup-ii].r ;
00119    }
00120    cf[nby2].r = 2.0 * row[nby2].r ; cf[nby2].i = 0.0 ;
00121    cg[nby2].r = 2.0 * row[nby2].i ; cg[nby2].i = 0.0 ;
00122 
00123    /* phase shift both rows (cf,cg) */
00124 
00125    dk = (2.0*PI) / nup ;
00126    sf = -af * dk ; sg = -ag * dk ;
00127 
00128 #ifdef RECUR
00129    csf = CEXPIT(sf) ; csg = CEXPIT(sg) ;
00130    fac.r = gac.r = 1.0 ;
00131    fac.i = gac.i = 0.0 ;
00132 #endif
00133 
00134    for( ii=1 ; ii <= nby2 ; ii++ ){
00135 #ifdef RECUR
00136       fac = CMULT( csf , fac ) ; cf[ii] = CMULT( fac , cf[ii] ) ;
00137       gac = CMULT( csg , gac ) ; cg[ii] = CMULT( gac , cg[ii] ) ;
00138 #else
00139       fac = CEXPIT(ii*sf) ; cf[ii] = CMULT( fac , cf[ii] ) ;
00140       gac = CEXPIT(ii*sg) ; cg[ii] = CMULT( gac , cg[ii] ) ;
00141 #endif
00142    }
00143    cf[nby2].i = 0.0 ; cg[nby2].i = 0.0 ;
00144 
00145    /* retangle the coefficients from 2 rows */
00146 
00147    row[0].r = cf[0].r ; row[0].i = cg[0].r ;
00148    for( ii=1 ; ii < nby2 ; ii++ ){
00149       row[ii].r     =  cf[ii].r - cg[ii].i ;
00150       row[ii].i     =  cf[ii].i + cg[ii].r ;
00151       row[nup-ii].r =  cf[ii].r + cg[ii].i ;
00152       row[nup-ii].i = -cf[ii].i + cg[ii].r ;
00153    }
00154    row[nby2].r = cf[nby2].r ;
00155    row[nby2].i = cg[nby2].r ;
00156 
00157    /* inverse FFT and store back in output arrays */
00158 
00159    csfft_cox( 1 , nup , row ) ;
00160 
00161    sf = 0.5 / nup ;              /* 0.5 to allow for twice too big above */
00162 
00163    if( g != NULL )
00164       for( ii=0; ii < n; ii++ ){ f[ii] = sf*row[ii].r; g[ii] = sf*row[ii].i; }
00165    else
00166       for( ii=0; ii < n; ii++ ){ f[ii] = sf*row[ii].r; }
00167 
00168    EXRETURN ;
00169 }
00170 
00171 /*--------------------------------------------------------------------------
00172   Some stuff needed for polynomial interpolation
00173 ----------------------------------------------------------------------------*/
00174 
00175 static int    nlcbuf = 0 ;     /* workspace */
00176 static float * lcbuf = NULL ;
00177 
00178    /* f[i], but inside the legal range */
00179 
00180 #ifdef ZFILL
00181 #  define FINS(i) ( ((i)<0 || (i)>=n) ? 0.0 : f[(i)] )
00182 #else
00183 #  define FINS(i) ( ((i)<0) ? f[0] : ((i)>=n) ? f[n-1] : f[(i)] )
00184 #endif
00185 
00186 #define SEPARATE_FINS
00187 
00188 /*---------------------------------------------------------------------------*/
00189 
00190 /*---------------------------------------------------------------------------
00191   Shift 1 row with with heptic Lagrange polynomial interpolation [Nov 1998].
00192   Note that heptic interpolation is about the same as Hamming-weighted
00193   3-sidelobe sinc interpolation .
00194 -----------------------------------------------------------------------------*/
00195 
00196    /* seventh order polynomials */
00197 
00198 #define S_M3(x) (x*(x*x-1.0)*(x*x-4.0)*(x-3.0)*(4.0-x)*0.0001984126984)
00199 #define S_M2(x) (x*(x*x-1.0)*(x-2.0)*(x*x-9.0)*(x-4.0)*0.001388888889)
00200 #define S_M1(x) (x*(x-1.0)*(x*x-4.0)*(x*x-9.0)*(4.0-x)*0.004166666667)
00201 #define S_00(x) ((x*x-1.0)*(x*x-4.0)*(x*x-9.0)*(x-4.0)*0.006944444444)
00202 #define S_P1(x) (x*(x+1.0)*(x*x-4.0)*(x*x-9.0)*(4.0-x)*0.006944444444)
00203 #define S_P2(x) (x*(x*x-1.0)*(x+2.0)*(x*x-9.0)*(x-4.0)*0.004166666667)
00204 #define S_P3(x) (x*(x*x-1.0)*(x*x-4.0)*(x+3.0)*(4.0-x)*0.001388888889)
00205 #define S_P4(x) (x*(x*x-1.0)*(x*x-4.0)*(x*x-9.0)*0.0001984126984)
00206 
00207 void hept_shift( int n , float af , float * f )
00208 {
00209    int   ii , ia , ix ;
00210    float  wt_m1,wt_00,wt_p1,wt_p2 , aa , wt_m2,wt_p3,wt_m3,wt_p4;
00211 #ifdef SEPARATE_FINS
00212    int ibot,itop ;
00213 #endif
00214 
00215 ENTRY("hept_shift") ;
00216 
00217    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00218 
00219    /* 15 Mar 2001: if shift is too large, return all zeros */
00220 
00221    if( ia <= -n || ia >= n ){
00222       for( ii=0 ; ii < n ; ii++ ) f[ii] = 0.0 ;
00223       EXRETURN ;
00224    }
00225 
00226    aa = af - ia ;
00227    wt_m1 = S_M1(aa) ; wt_00 = S_00(aa) ;
00228    wt_p1 = S_P1(aa) ; wt_p2 = S_P2(aa) ;
00229    wt_m2 = S_M2(aa) ; wt_p3 = S_P3(aa) ;
00230    wt_m3 = S_M3(aa) ; wt_p4 = S_P4(aa) ;
00231 
00232    if( n > nlcbuf ){
00233       if( lcbuf != NULL ) free(lcbuf) ;
00234       lcbuf  = (float *) malloc( sizeof(float) * n ) ;
00235       nlcbuf = n ;
00236    }
00237 
00238 #ifdef SEPARATE_FINS
00239    ibot = 3-ia ;   if( ibot < 0   ) ibot = 0 ;
00240    itop = n-5-ia ; if( itop > n-1 ) itop = n-1 ;
00241 
00242    for( ii=ibot ; ii <= itop ; ii++ ){
00243       ix = ii + ia ;
00244       lcbuf[ii] =  wt_m2 * f[ix-2] + wt_m1 * f[ix-1] + wt_00 * f[ix]
00245                  + wt_p1 * f[ix+1] + wt_p2 * f[ix+2] + wt_p3 * f[ix+3]
00246                  + wt_m3 * f[ix-3] + wt_p4 * f[ix+4] ;
00247    }
00248 
00249    if( ibot > n ) ibot = n ; /* 15 Mar 2001 */
00250    for( ii=0 ; ii < ibot ; ii++ ){
00251       ix = ii + ia ;
00252       lcbuf[ii] =  wt_m2 * FINS(ix-2) + wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00253                  + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) + wt_p3 * FINS(ix+3)
00254                  + wt_m3 * FINS(ix-3) + wt_p4 * FINS(ix+4) ;
00255    }
00256 
00257    if( itop < 0 ) itop = -1 ; /* 15 Mar 2001 */
00258    for( ii=itop+1 ; ii < n ; ii++ ){
00259       ix = ii + ia ;
00260       lcbuf[ii] =  wt_m2 * FINS(ix-2) + wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00261                  + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) + wt_p3 * FINS(ix+3)
00262                  + wt_m3 * FINS(ix-3) + wt_p4 * FINS(ix+4) ;
00263    }
00264 #else /* not SEPARATE_FINS */
00265    for( ii=0 ; ii < n ; ii++ ){
00266       ix = ii + ia ;
00267       if( ix > 1 && ix < n-3 )
00268          lcbuf[ii] =  wt_m2 * f[ix-2] + wt_m1 * f[ix-1] + wt_00 * f[ix]
00269                     + wt_p1 * f[ix+1] + wt_p2 * f[ix+2] + wt_p3 * f[ix+3]
00270                     + wt_m3 * f[ix-3] + wt_p4 * f[ix+4] ;
00271       else
00272          lcbuf[ii] =  wt_m2 * FINS(ix-2) + wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00273                     + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) + wt_p3 * FINS(ix+3)
00274                     + wt_m3 * FINS(ix-3) + wt_p4 * FINS(ix+4) ;
00275    }
00276 #endif /* SEPARATE_FINS */
00277 
00278    memcpy( f , lcbuf , sizeof(float)*n ) ;
00279    EXRETURN ;
00280 }
00281 
00282 void hept_shift2( int n, int nup, float af, float * f, float ag, float * g )
00283 {
00284                    hept_shift( n , af , f ) ;
00285    if( g != NULL ) hept_shift( n , ag , g ) ;
00286    return ;
00287 }
00288 
00289 /*---------------------------------------------------------------------------
00290   Shift 1 row with with quintic Lagrange polynomial interpolation [Nov 1998].
00291   Note that quintic interpolation is about the same as Hamming-weighted
00292   2-sidelobe sinc interpolation .
00293 -----------------------------------------------------------------------------*/
00294 
00295    /* quintic interpolation polynomials (Lagrange) */
00296 
00297 #define Q_M2(x)  (x*(x*x-1.0)*(2.0-x)*(x-3.0)*0.008333333)
00298 #define Q_M1(x)  (x*(x*x-4.0)*(x-1.0)*(x-3.0)*0.041666667)
00299 #define Q_00(x)  ((x*x-4.0)*(x*x-1.0)*(3.0-x)*0.083333333)
00300 #define Q_P1(x)  (x*(x*x-4.0)*(x+1.0)*(x-3.0)*0.083333333)
00301 #define Q_P2(x)  (x*(x*x-1.0)*(x+2.0)*(3.0-x)*0.041666667)
00302 #define Q_P3(x)  (x*(x*x-1.0)*(x*x-4.0)*0.008333333)
00303 
00304 void quint_shift( int n , float af , float * f )
00305 {
00306    int   ii , ia , ix ;
00307    float  wt_m1 , wt_00 , wt_p1 , wt_p2 , aa , wt_m2 , wt_p3 ;
00308 #ifdef SEPARATE_FINS
00309    int ibot,itop ;
00310 #endif
00311 
00312 ENTRY("quint_shift") ;
00313 
00314    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00315 
00316    /* 15 Mar 2001: if shift is too large, return all zeros */
00317 
00318    if( ia <= -n || ia >= n ){
00319       for( ii=0 ; ii < n ; ii++ ) f[ii] = 0.0 ;
00320       EXRETURN ;
00321    }
00322 
00323    aa = af - ia ;
00324    wt_m1 = Q_M1(aa) ; wt_00 = Q_00(aa) ;
00325    wt_p1 = Q_P1(aa) ; wt_p2 = Q_P2(aa) ;
00326    wt_m2 = Q_M2(aa) ; wt_p3 = Q_P3(aa) ;
00327 
00328    if( n > nlcbuf ){
00329       if( lcbuf != NULL ) free(lcbuf) ;
00330       lcbuf  = (float *) malloc( sizeof(float) * n ) ;
00331       nlcbuf = n ;
00332    }
00333 
00334 #ifdef SEPARATE_FINS
00335    ibot = 2-ia ;   if( ibot < 0   ) ibot = 0 ;
00336    itop = n-4-ia ; if( itop > n-1 ) itop = n-1 ;
00337 
00338    for( ii=ibot ; ii <= itop ; ii++ ){
00339       ix = ii + ia ;
00340       lcbuf[ii] =  wt_m2 * f[ix-2] + wt_m1 * f[ix-1] + wt_00 * f[ix]
00341                  + wt_p1 * f[ix+1] + wt_p2 * f[ix+2] + wt_p3 * f[ix+3] ;
00342    }
00343 
00344    if( ibot > n ) ibot = n ; /* 15 Mar 2001 */
00345    for( ii=0 ; ii < ibot ; ii++ ){
00346       ix = ii + ia ;
00347       lcbuf[ii] =  wt_m2 * FINS(ix-2) + wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00348                  + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) + wt_p3 * FINS(ix+3) ;
00349    }
00350 
00351    if( itop < 0 ) itop = -1 ; /* 15 Mar 2001 */
00352    for( ii=itop+1 ; ii < n ; ii++ ){
00353       ix = ii + ia ;
00354       lcbuf[ii] =  wt_m2 * FINS(ix-2) + wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00355                  + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) + wt_p3 * FINS(ix+3) ;
00356    }
00357 #else /* not SEPARATE_FINS */
00358    for( ii=0 ; ii < n ; ii++ ){
00359       ix = ii + ia ;
00360       if( ix > 1 && ix < n-3 )
00361          lcbuf[ii] =  wt_m2 * f[ix-2] + wt_m1 * f[ix-1] + wt_00 * f[ix]
00362                     + wt_p1 * f[ix+1] + wt_p2 * f[ix+2] + wt_p3 * f[ix+3] ;
00363       else
00364          lcbuf[ii] =  wt_m2 * FINS(ix-2) + wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00365                     + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) + wt_p3 * FINS(ix+3) ;
00366    }
00367 #endif /* SEPARATE_FINS */
00368 
00369    memcpy( f , lcbuf , sizeof(float)*n ) ;
00370    EXRETURN ;
00371 }
00372 
00373 void quint_shift2( int n, int nup, float af, float * f, float ag, float * g )
00374 {
00375                    quint_shift( n , af , f ) ;
00376    if( g != NULL ) quint_shift( n , ag , g ) ;
00377    return ;
00378 }
00379 
00380 /*---------------------------------------------------------------------------
00381   Shift 1 row with with cubic interpolation
00382 -----------------------------------------------------------------------------*/
00383 
00384    /* cubic interpolation polynomials */
00385 
00386 #define P_M1(x)  ((x)*(1.0-(x))*((x)-2.0)*0.1666667)
00387 #define P_00(x)  (((x)+1.0)*((x)-1.0)*((x)-2.0)*0.5)
00388 #define P_P1(x)  ((x)*((x)+1.0)*(2.0-(x))*0.5)
00389 #define P_P2(x)  ((x)*((x)+1.0)*((x)-1.0)*0.1666667)
00390 
00391 void cub_shift( int n , float af , float * f )
00392 {
00393    int   ii , ia , ix ;
00394    float  wt_m1 , wt_00 , wt_p1 , wt_p2 , aa ;
00395 #ifdef SEPARATE_FINS
00396    int ibot,itop ;
00397 #endif
00398 
00399 ENTRY("cub_shift") ;
00400 
00401    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00402 
00403    /* 15 Mar 2001: if shift is too large, return all zeros */
00404 
00405    if( ia <= -n || ia >= n ){
00406       for( ii=0 ; ii < n ; ii++ ) f[ii] = 0.0 ;
00407       EXRETURN ;
00408    }
00409 
00410    aa = af - ia ;
00411    wt_m1 = P_M1(aa) ; wt_00 = P_00(aa) ;
00412    wt_p1 = P_P1(aa) ; wt_p2 = P_P2(aa) ;
00413 
00414    if( n > nlcbuf ){
00415       if( lcbuf != NULL ) free(lcbuf) ;
00416       lcbuf  = (float *) malloc( sizeof(float) * n ) ;
00417       nlcbuf = n ;
00418    }
00419 
00420 #ifdef SEPARATE_FINS
00421    ibot = 1-ia ;   if( ibot < 0   ) ibot = 0 ;
00422    itop = n-3-ia ; if( itop > n-1 ) itop = n-1 ;
00423 
00424    for( ii=ibot ; ii <= itop ; ii++ ){
00425       ix = ii + ia ;
00426       lcbuf[ii] =  wt_m1 * f[ix-1] + wt_00 * f[ix]
00427                  + wt_p1 * f[ix+1] + wt_p2 * f[ix+2] ;
00428    }
00429 
00430    if( ibot > n ) ibot = n ; /* 15 Mar 2001 */
00431    for( ii=0 ; ii < ibot ; ii++ ){
00432       ix = ii + ia ;
00433       lcbuf[ii] =  wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00434                  + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) ;
00435    }
00436 
00437    if( itop < 0 ) itop = -1 ; /* 15 Mar 2001 */
00438    for( ii=itop+1 ; ii < n ; ii++ ){
00439       ix = ii + ia ;
00440       lcbuf[ii] =  wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00441                  + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) ;
00442    }
00443 #else /* not SEPARATE_FINS */
00444    for( ii=0 ; ii < n ; ii++ ){
00445       ix = ii + ia ;
00446       if( ix > 0 && ix < n-2 )
00447          lcbuf[ii] =  wt_m1 * f[ix-1] + wt_00 * f[ix]
00448                     + wt_p1 * f[ix+1] + wt_p2 * f[ix+2] ;
00449       else
00450          lcbuf[ii] =  wt_m1 * FINS(ix-1) + wt_00 * FINS(ix)
00451                     + wt_p1 * FINS(ix+1) + wt_p2 * FINS(ix+2) ;
00452    }
00453 #endif /* SEPARATE_FINS */
00454 
00455    memcpy( f , lcbuf , sizeof(float)*n ) ;
00456    EXRETURN ;
00457 }
00458 
00459 void cub_shift2( int n, int nup, float af, float * f, float ag, float * g )
00460 {
00461                    cub_shift( n , af , f ) ;
00462    if( g != NULL ) cub_shift( n , ag , g ) ;
00463    return ;
00464 }
00465 
00466 /*---------------------------------------------------------------------------
00467    Shift one row with linear interpolation
00468 -----------------------------------------------------------------------------*/
00469 
00470 void lin_shift( int n , float af , float * f )
00471 {
00472    int   ii , ia , ix ;
00473    float  wt_00 , wt_p1 , aa ;
00474 #ifdef SEPARATE_FINS
00475    int ibot,itop ;
00476 #endif
00477 
00478 ENTRY("lin_shift") ;
00479 
00480    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00481    aa = af - ia ;
00482    wt_00 = 1.0 - aa ; wt_p1 = aa ;  /* linear interpolation weights */
00483 
00484    /* 15 Mar 2001: if shift is too large, return all zeros */
00485 
00486    if( ia <= -n || ia >= n ){
00487       for( ii=0 ; ii < n ; ii++ ) f[ii] = 0.0 ;
00488       EXRETURN ;
00489    }
00490 
00491    if( n > nlcbuf ){
00492       if( lcbuf != NULL ) free(lcbuf) ;
00493       lcbuf  = (float *) malloc( sizeof(float) * n ) ;
00494       nlcbuf = n ;
00495    }
00496 
00497 #ifdef SEPARATE_FINS
00498    ibot = -ia  ;   if( ibot < 0   ) ibot = 0 ;
00499    itop = n-2-ia ; if( itop > n-1 ) itop = n-1 ;
00500 
00501 #if 0
00502 if(PRINT_TRACING){
00503   char str[256]; sprintf(str,"n=%d ia=%d ibot=%d itop=%d",n,ia,ibot,itop); STATUS(str);
00504 }
00505 #endif
00506 
00507    for( ii=ibot ; ii <= itop ; ii++ ){
00508       ix = ii + ia ;
00509       lcbuf[ii] =  wt_00 * f[ix] + wt_p1 * f[ix+1] ;
00510    }
00511 
00512    if( ibot > n ) ibot = n ; /* 15 Mar 2001 */
00513    for( ii=0 ; ii < ibot ; ii++ ){
00514       ix = ii + ia ;
00515       lcbuf[ii] =  wt_00 * FINS(ix) + wt_p1 * FINS(ix+1) ;
00516    }
00517 
00518    if( itop < 0 ) itop = -1 ; /* 15 Mar 2001 */
00519    for( ii=itop+1 ; ii < n ; ii++ ){
00520       ix = ii + ia ;
00521       lcbuf[ii] =  wt_00 * FINS(ix) + wt_p1 * FINS(ix+1) ;
00522    }
00523 #else
00524    for( ii=0 ; ii < n ; ii++ ){
00525       ix = ii + ia ;
00526       if( ix >= 0 && ix < n-1 )
00527          lcbuf[ii] =  wt_00 * f[ix] + wt_p1 * f[ix+1] ;
00528       else
00529          lcbuf[ii] =  wt_00 * FINS(ix) + wt_p1 * FINS(ix+1) ;
00530    }
00531 #endif /* SEPARATE_FINS */
00532 
00533    memcpy( f , lcbuf , sizeof(float)*n ) ;
00534    EXRETURN ;
00535 }
00536 
00537 void lin_shift2( int n, int nup, float af, float * f, float ag, float * g )
00538 {
00539                    lin_shift( n , af , f ) ;
00540    if( g != NULL ) lin_shift( n , ag , g ) ;
00541    return ;
00542 }
00543 
00544 /*--------------------------------------------------------------------------
00545   This is for experimental purposes only -- DO NOT USE!
00546 ----------------------------------------------------------------------------*/
00547 
00548 void nn_shift( int n , float af , float * f )
00549 {
00550    int   ii , ia , ix ;
00551 
00552 ENTRY("nn_shift") ;
00553 
00554    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00555 
00556    /* 15 Mar 2001: if shift is too large, return all zeros */
00557 
00558    if( ia <= -n || ia >= n ){
00559       for( ii=0 ; ii < n ; ii++ ) f[ii] = 0.0 ;
00560       EXRETURN ;
00561    }
00562 
00563    if( n > nlcbuf ){
00564       if( lcbuf != NULL ) free(lcbuf) ;
00565       lcbuf  = (float *) malloc( sizeof(float) * n ) ;
00566       nlcbuf = n ;
00567    }
00568 
00569    for( ii=0 ; ii < n ; ii++ ){
00570       ix = ii + ia ;
00571       lcbuf[ii] = FINS(ix) ;
00572    }
00573 
00574    memcpy( f , lcbuf , sizeof(float)*n ) ;
00575    EXRETURN ;
00576 }
00577 
00578 void nn_shift2( int n, int nup, float af, float * f, float ag, float * g )
00579 {
00580                    nn_shift( n , af , f ) ;
00581    if( g != NULL ) nn_shift( n , ag , g ) ;
00582    return ;
00583 }
00584 
00585 /*---------------------------------------------------------------------------
00586    More experiments: two-step interpolation
00587 -----------------------------------------------------------------------------*/
00588 
00589 void ts_shift( int n , float af , float * f )
00590 {
00591    register int ii , ia , ix ;
00592    float aa ;
00593    int ibot,itop ;
00594 
00595    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00596 
00597    /* 15 Mar 2001: if shift is too large, return all zeros */
00598 
00599    if( ia <= -n || ia >= n ){
00600       for( ii=0 ; ii < n ; ii++ ) f[ii] = 0.0 ;
00601       EXRETURN ;
00602    }
00603 
00604    aa = af - ia ;
00605 
00606    if( n > nlcbuf ){
00607       if( lcbuf != NULL ) free(lcbuf) ;
00608       lcbuf  = (float *) malloc( sizeof(float) * n ) ;
00609       nlcbuf = n ;
00610    }
00611 
00612    ibot = -ia  ;   if( ibot < 0   ) ibot = 0 ;
00613    itop = n-2-ia ; if( itop > n-1 ) itop = n-1 ;
00614 
00615    if( aa < 0.30 ){
00616       memcpy( lcbuf+ibot, f+(ibot+ia)  , (itop+1-ibot)*sizeof(float) );
00617       for( ii=0 ; ii < ibot ; ii++ ){
00618          ix = ii + ia ; lcbuf[ii] = FINS(ix) ;
00619       }
00620       for( ii=itop+1 ; ii < n ; ii++ ){
00621          ix = ii + ia ; lcbuf[ii] = FINS(ix) ;
00622       }
00623    }
00624 
00625    else if( aa > 0.70 ){
00626       memcpy( lcbuf+ibot, f+(ibot+1+ia), (itop+1-ibot)*sizeof(float) );
00627       for( ii=0 ; ii < ibot ; ii++ ){
00628          ix = ii + ia ; lcbuf[ii] = FINS(ix+1) ;
00629       }
00630       for( ii=itop+1 ; ii < n ; ii++ ){
00631          ix = ii + ia ; lcbuf[ii] = FINS(ix+1) ;
00632       }
00633 
00634    } else {
00635       for( ii=ibot ; ii <= itop ; ii++ ){
00636          ix = ii + ia ; lcbuf[ii] =  0.5*( f[ix] + f[ix+1] ) ;
00637       }
00638       if( ibot > n ) ibot = n ; /* 15 Mar 2001 */
00639       for( ii=0 ; ii < ibot ; ii++ ){
00640          ix = ii + ia ; lcbuf[ii] =  0.5*( FINS(ix) + FINS(ix+1) ) ;
00641       }
00642       if( itop < 0 ) itop = -1 ; /* 15 Mar 2001 */
00643       for( ii=itop+1 ; ii < n ; ii++ ){
00644          ix = ii + ia ; lcbuf[ii] =  0.5*( FINS(ix) + FINS(ix+1) ) ;
00645       }
00646    }
00647    memcpy( f , lcbuf , sizeof(float)*n ) ;
00648    return ;
00649 }
00650 
00651 void ts_shift2( int n, int nup, float af, float * f, float ag, float * g )
00652 {
00653                    ts_shift( n , af , f ) ;
00654    if( g != NULL ) ts_shift( n , ag , g ) ;
00655    return ;
00656 }