/** * Warps a sampled function using only 1st components of shifts. * @param f array of values f(x). * @return array of values g(y) = f(y-u1(x(y)). */ public float[][] warp1(float[][] f) { SincInterpolator si = new SincInterpolator(); float[][] g = new float[_n2][_n1]; for (int i2=0; i2<_n2; ++i2) { for (int i1=0; i1<_n1; ++i1) { g[i2][i1] = si.interpolate(_n1,1.0,0.0,f[i2],i1-u1y(i1,i2)); } } return g; }
/** * Accumulates a specified real value y(x) into uniformly-sampled yu. * Accumulation is the transpose (not the inverse) of interpolation. * Whereas interpolation gathers from uniformly sampled values. * accumulation scatters into uniformly sampled values. * @param xa value x at which to accumulate. * @param ya value y(x) to accumulate. * @param nxu number of input/output samples. * @param dxu input/output sampling interval. * @param fxu first input/output sampled x value. * @param yu input/output array of sampled values y(x). */ public void accumulate( double xa, float ya, int nxu, double dxu, double fxu, float[] yu) { double xscale = 1.0/dxu; double xshift = _lsinc-fxu*xscale; int nxum = nxu-_lsinc; accumulate(xscale,xshift,nxum,xa,ya,nxu,yu); }
/** * Interpolates multiple complex values y(x). * Complex output samples are packed in the specified output array as * real, imag, real, imag, .... * @param nxu number of input samples. * @param dxu input sampling interval. * @param fxu first input sampled x value. * @param yu input array[2*nxu] of sampled complex values y(x). * @param nxi number of output samples. * @param xi input array[nxi] of x values at which to interpolate. * @param yi output array[2*nxi] of interpolated complex values y(x). */ public void interpolateComplex( int nxu, double dxu, double fxu, float[] yu, int nxi, float[] xi, float[] yi) { double xscale = 1.0/dxu; double xshift = _lsinc-fxu*xscale; int nxum = nxu-_lsinc; for (int ixi=0; ixi<nxi; ++ixi) interpolateComplex(xscale,xshift,nxum,nxu,yu,ixi,xi[ixi],yi); }
/** * Interpolates multiple real values y(x). * @param nxu number of input samples. * @param dxu input sampling interval. * @param fxu first input sampled x value. * @param yu input array of sampled values y(x). * @param nxi number of output samples. * @param dxi output sampling interval. * @param fxi first output sampled x value. * @param yi output array of interpolated values y(x). */ public void interpolate( int nxu, double dxu, double fxu, float[] yu, int nxi, double dxi, double fxi, float[] yi) { if (dxu==dxi) { shift(nxu,dxu,fxu,yu,nxi,fxi,yi); } else { double xscale = 1.0/dxu; double xshift = _lsinc-fxu*xscale; int nxum = nxu-_lsinc; for (int ixi=0; ixi<nxi; ++ixi) { double xi = fxi+ixi*dxi; yi[ixi] = interpolate(xscale,xshift,nxum,nxu,yu,xi); } } }
public void testComplex() { SincInterpolator si = new SincInterpolator(); Random random = new Random(); int nxu = 100; float[] zi = new float[nx]; float[] zc = new float[2*nx]; si.setExtrapolation(SincInterpolator.Extrapolation.ZERO); si.interpolate(nxu,dxu,fxu,yr,nx,dx,fx,zr); si.interpolate(nxu,dxu,fxu,yi,nx,dx,fx,zi); si.interpolateComplex(nxu,dxu,fxu,yc,nx,dx,fx,zc); for (int ix=0; ix<nx; ++ix) { assertEquals(zr[ix],zc[2*ix ],0.0); assertEquals(zi[ix],zc[2*ix+1],0.0); si.setExtrapolation(SincInterpolator.Extrapolation.CONSTANT); si.interpolate(nxu,dxu,fxu,yr,nx,dx,fx,zr); si.interpolate(nxu,dxu,fxu,yi,nx,dx,fx,zi); si.interpolateComplex(nxu,dxu,fxu,yc,nx,dx,fx,zc); for (int ix=0; ix<nx; ++ix) { assertEquals(zr[ix],zc[2*ix ],0.0);
SincInterpolator si = new SincInterpolator(); si.setExtrapolation(extrapolation); si.interpolate(nxu, dxu, fxu, yu, nx, x, yi); si.accumulate(nx, x, y, nxu, dxu, fxu, ya);
private void testInterpolatorWithSweep(SincInterpolator si) { double emax = si.getMaximumError(); double fmax = si.getMaximumFrequency(); int lmax = si.getMaximumLength(); long nbytes = si.getTableBytes(); trace("lmax="+lmax+" fmax="+fmax+" emax="+emax+" nbytes="+nbytes); yu[ixu] = (float)sweep(fmax,nmax,x); si.setExtrapolation(SincInterpolator.Extrapolation.CONSTANT); dx = (xmax-fx)/(nx-1); float[] y = new float[nx]; si.interpolate(nxu,dxu,fxu,yu,nx,dx,fx,y); dx = dxu; fx = fxu+shift; si.interpolate(nxu,dxu,fxu,yu,nx,dx,fx,y); error = 0.0; for (int ix=0; ix<nx; ++ix) {
/** * Computes a sequence warped by applying specified shifts. * @param u input array of shifts. * @param g input array for the sequence to be warped. * @param h output array for the warped sequence. */ public void applyShifts(float[] u, float[] g, float[] h) { int n1 = u.length; for (int i1=0; i1<n1; ++i1) { h[i1] = _si.interpolate(n1,1.0,0.0,g,i1+u[i1]); } }
public void testExtrapolation() { SincInterpolator si = new SincInterpolator(); Random random = new Random(); int lmax = si.getMaximumLength(); int nxu = 2*lmax; int npad = lmax; yc[npad+nxu+ipad] = yc[npad+nxu-1]; si.setExtrapolation(SincInterpolator.Extrapolation.ZERO); si.interpolate(nxu,dxu,fxu,yi,nx,dx,fx,yo); si.interpolate(npad+nxu+npad,dxu,0.0,yz,nx,dx,fx,yt); for (int ix=0; ix<nx; ++ix) assertEquals(yo[ix],yt[ix],0.0); si.setExtrapolation(SincInterpolator.Extrapolation.CONSTANT); si.interpolate(nxu,dxu,fxu,yi,nx,dx,fx,yo); si.interpolate(npad+nxu+npad,dxu,0.0,yc,nx,dx,fx,yt); for (int ix=0; ix<nx; ++ix) assertEquals(yo[ix],yt[ix],0.0);
SincInterpolator si = new SincInterpolator(); si.setExtrapolation(SincInterpolator.Extrapolation.CONSTANT); for (int i2=0; i2<nl2; ++i2) { for (int i1=0; i1<nl1; ++i1) { spyr[lev][0][i2][i1] += si.interpolate( m1,2,0,m2,2,0,spyr[lev+1][0],i1,i2);
float[] lo12 = zerofloat(m2); float[] lo13 = zerofloat(m3); SincInterpolator si = SincInterpolator.fromErrorAndFrequency(0.001,0.4); si.setExtrapolation(SincInterpolator.Extrapolation.CONSTANT); for (int i3=0; i3<nl3; i3=i3+2) { j3 = i3/2; spyr[lev][1][i3][i2][i1] = si.interpolate(m1,2,0,lo11,i1); spyr[lev][1][i3][i2][i1] = si.interpolate(m2,2,0,lo12,i2); spyr[lev][1][i3][i2][i1] = si.interpolate(m3,2,0,lo13,i3);
/** * Construct a shift estimator with specified parameters. * When applied to multi-dimensional arrays, the estimator has half-width * sigma1 for the 1st dimension, half-width sigma2 for the 2nd dimension, * and half-width sigma3 for 3rd and higher dimensions. * @param sigma1 correlation window half-width for 1st dimension; * must not be less than 1. * @param sigma2 correlation window half-width for 2nd dimension; * must not be less than 1. * @param sigma3 correlation window half-width for 3rd and higher * dimensions; must not be less than 1. */ public LocalShiftFinder(double sigma1, double sigma2, double sigma3) { Check.argument(sigma1>=1.0,"sigma1>=1.0"); Check.argument(sigma2>=1.0,"sigma2>=1.0"); Check.argument(sigma3>=1.0,"sigma3>=1.0"); _lcfSimple = new LocalCorrelationFilter( LocalCorrelationFilter.Type.SIMPLE, LocalCorrelationFilter.Window.GAUSSIAN, sigma1,sigma2,sigma3); _si = new SincInterpolator(); _si.setExtrapolation(SincInterpolator.Extrapolation.CONSTANT); }
/** * Returns a sinc interpolator with specified maximum error and frequency. * Computes the maximum length lmax. * @param emax the maximum error for frequencies less than fmax; e.g., * 0.01 for 1% percent error. Must be greater than 0.0 and less than 1.0. * @param fmax the maximum frequency, in cycles per sample. * Must be greater than 0.0 and less than 0.5. * @return the sinc interpolator. */ public static SincInterpolator fromErrorAndFrequency( double emax, double fmax) { return new SincInterpolator(emax,fmax,0); }
public void testErrorAndFrequency() { for (double emax:_emaxs) { for (double fmax:_fmaxs) { SincInterpolator si = SincInterpolator.fromErrorAndFrequency(emax,fmax); testInterpolator(si); } } }
public void testFrequencyAndLength() { for (double fmax:_fmaxs) { for (int lmax:_lmaxs) { if ((1.0-2.0*fmax)*lmax>1.0) { SincInterpolator si = SincInterpolator.fromFrequencyAndLength(fmax,lmax); testInterpolator(si); } } } }
public void testErrorAndLength() { for (double emax:_emaxs) { for (int lmax:_lmaxs) { SincInterpolator si = SincInterpolator.fromErrorAndLength(emax,lmax); testInterpolator(si); } } }
public void compute(int i2) { for (int i1=0; i1<n1; ++i1) { hf[i2][i1] = _si.interpolate(n1,1.0,0.0,gf[i2],i1+uf[i2][i1]); } }}); }
/** * Returns a sinc interpolator with specified maximum error and length. * Computes the maximum frequency fmax. Note that for some parameters * emax and lmax, the maximum frequency fmax may be zero. In this case, * the returned interpolator is useless. * @param emax the maximum error for frequencies less than fmax; e.g., * 0.01 for 1% percent error. 0.0 < emax <= 0.1 is required. * @param lmax the maximum interpolator length, in samples. * Must be an even integer not less than 8. * @return the sinc interpolator. */ public static SincInterpolator fromErrorAndLength( double emax, int lmax) { return new SincInterpolator(emax,0.0,lmax); }
/** * Unwarps a sampled function using only 1st components of shifts. * @param g array of values g(x). * @return array of values f(x) = g(x+u1(x)). */ public float[][] unwarp1(float[][] g) { SincInterpolator si = new SincInterpolator(); float[][] f = new float[_n2][_n1]; for (int i2=0; i2<_n2; ++i2) { for (int i1=0; i1<_n1; ++i1) { f[i2][i1] = si.interpolate(_n1,1.0,0.0,g[i2],i1+u1x(i1,i2)); } } return f; }
/** * Interpolates one real value y(x). * @param nxu number of input samples. * @param dxu input sampling interval. * @param fxu first input sampled x value. * @param yu input array of sampled values y(x). * @param xi value x at which to interpolate. * @return interpolated value y(x). */ public float interpolate( int nxu, double dxu, double fxu, float[] yu, double xi) { double xscale = 1.0/dxu; double xshift = _lsinc-fxu*xscale; int nxum = nxu-_lsinc; return interpolate(xscale,xshift,nxum,nxu,yu,xi); }