edu.mines.jtk.dsp.SincInterp Maven / Gradle / Ivy
Show all versions of edu-mines-jtk Show documentation
/****************************************************************************
Copyright 2012, Colorado School of Mines and others.
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
****************************************************************************/
package edu.mines.jtk.dsp;
import java.util.HashMap;
import static edu.mines.jtk.util.ArrayMath.*;
import edu.mines.jtk.util.Check;
/**
* A sinc interpolator for bandlimited uniformly-sampled functions y(x).
* Interpolators can be designed for any two of three parameters: maximum
* error (emax), maximum frequency (fmax) and maximum length (lmax). The
* parameter not specified is computed when an interpolator is designed.
*
* Below the specified (or computed) maximum frequency fmax, the maximum
* interpolation error should be less than the specified (or computed)
* maximum error emax. For frequencies above fmax, interpolation error
* may be much greater. Therefore, sequences to be interpolated should
* be bandlimited to frequencies less than fmax.
*
* The maximum length lmax of an interpolator is an even positive integer.
* It is the number of uniform samples required to interpolate a single
* value y(x). Ideally, the weights applied to each uniform sample are
* values of a sinc function. Although the ideal sinc function yields zero
* interpolation error for all frequencies up to the Nyquist frequency
* (0.5 cycles/sample), it has infinite length.
*
* With recursive filtering, infinite-length approximations to the sinc
* function are feasible and, in some applications, most efficient. When
* the number of interpolated values is large relative to the number of
* uniform samples, the cost of recursive filtering is amortized over those
* many interpolated values, and can be negligible. However, this cost
* becomes significant when only a few values are interpolated for each
* sequence of uniform samples.
*
* This interpolator is based on a finite-length approximation
* to the sinc function. The efficiency of finite-length interpolators
* like this one does not depend on the number of samples interpolated.
* Also, this interpolator is robust in the presence of noise spikes,
* which affect only nearby samples.
*
* Finite-length interpolators present a tradeoff between cost and accuracy.
* Interpolators with small maximum lengths are most efficient, and those
* with high maximum frequencies and small maximum errors are most accurate.
*
* When interpolating multiple values of y(x) from a single sequence of
* uniformly sampled values, efficiency may be improved by using one of the
* methods that enables specification of multiple x values at which to
* interpolate.
*
* @author Dave Hale, Colorado School of Mines
* @author Bill Harlan, Landmark Graphics
* @version 2012.12.21
* @deprecated Use class {@link edu.mines.jtk.dsp.SincInterpolator} instead.
*/
public class SincInterp {
/**
* The method used to extrapolate samples when interpolating.
* Sampled functions are defined implicitly by extrapolation outside the
* domain for which samples are specified explicitly, with either zero or
* constant values. If constant, the extrapolated values are the first and
* last uniform sample values. The default is extrapolation with zeros.
*/
public enum Extrapolation {
ZERO,
CONSTANT,
}
/**
* 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 SincInterp fromErrorAndLength(
double emax, int lmax)
{
return new SincInterp(emax,0.0,lmax);
}
/**
* 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 SincInterp fromErrorAndFrequency(
double emax, double fmax)
{
return new SincInterp(emax,fmax,0);
}
/**
* Returns a sinc interpolator with specified maximum frequency and length.
* Computes the maximum error emax.
*
* The product (1-2*fmax)*lmax must be greater than one. For when this
* product is less than one, a useful upper bound on interpolation error
* cannot be computed.
* @param fmax the maximum frequency, in cycles per sample.
* Must be greater than 0.0 and less than 0.5*(1.0-1.0/lmax).
* @param lmax the maximum interpolator length, in samples.
* Must be an even integer not less than 8 and greater than
* 1.0/(1.0-2.0*fmax).
* @return the sinc interpolator.
*/
public static SincInterp fromFrequencyAndLength(
double fmax, int lmax)
{
return new SincInterp(0.0,fmax,lmax);
}
/**
* Constructs a default sinc interpolator. The default design parameters
* are fmax = 0.3 cycles/sample (60% of Nyquist) and lmax = 8 samples.
* For these parameters, the computed maximum error is less than 0.007
* (0.7%). In testing, observed maximum error is less than 0.005 (0.5%).
*/
public SincInterp() {
this(0.0,0.3,8);
}
/**
* Gets the maximum error for this interpolator.
* @return the maximum error.
*/
public double getMaximumError() {
return _table.design.emax;
}
/**
* Gets the maximum frequency for this interpolator.
* @return the maximum frequency.
*/
public double getMaximumFrequency() {
return _table.design.fmax;
}
/**
* Gets the maximum length for this interpolator.
* @return the maximum length.
*/
public int getMaximumLength() {
return _table.design.lmax;
}
/**
* Gets the number of bytes consumed by the table of interpolators.
* The use of interpolators with small emax and large lmax may require
* the computation of large tables. This method can be used to determine
* how much memory is consumed by the table for an interpolator.
* @return the number of bytes.
*/
public long getTableBytes() {
long nbytes = 4L;
nbytes *= _table.lsinc;
nbytes *= _table.nsinc;
return nbytes;
}
/**
* Gets the extrapolation method for this interpolator.
* @return the extrapolation method.
*/
public Extrapolation getExtrapolation() {
return _extrap;
}
/**
* Sets the extrapolation method for this interpolator.
* The default extrapolation method is extrapolation with zeros.
* @param extrap the extrapolation method.
*/
public void setExtrapolation(Extrapolation extrap) {
_extrap = extrap;
}
/**
* 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);
}
/**
* 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 xi input array of x values at which to interpolate.
* @param yi output array of interpolated values y(x).
*/
public void interpolate(
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=8,"lmax>=8");
Check.argument(lmax%2==0,"lmax is even");
Check.argument((1.0-2.0*fmax)*lmax>1.0,"(1.0-2.0*fmax)*lmax>1.0");
} else if (fmax==0.0) {
Check.argument(emax<=0.1,"emax<=0.1");
Check.argument(lmax>=8,"lmax>=8");
Check.argument(lmax%2==0,"lmax is even");
} else if (lmax==0) {
Check.argument(emax<=0.1,"emax<=0.1");
Check.argument(fmax<0.5,"fmax<0.5");
}
_table = getTable(emax,fmax,lmax);
_lsinc = _table.lsinc;
_nsinc = _table.nsinc;
_nsincm1 = _table.nsincm1;
_ishift = _table.ishift;
_dsinc = _table.dsinc;
_asinc = _table.asinc;
}
/**
* Design parameters.
*/
private static class Design {
double emax;
double fmax;
int lmax;
Design(double emax, double fmax, int lmax) {
this.emax = emax;
this.fmax = fmax;
this.lmax = lmax;
}
public int hashCode() {
long lemax = Double.doubleToLongBits(emax);
long lfmax = Double.doubleToLongBits(fmax);
return (int)(lemax^(lemax>>>32)) ^
(int)(lfmax^(lfmax>>>32)) ^
lmax;
}
public boolean equals(Object object) {
Design that = (Design)object;
return this.emax==that.emax &&
this.fmax==that.fmax &&
this.lmax==that.lmax;
}
}
/**
* Table of sinc interpolator coefficients.
*/
private static class Table {
Design design; // here, all three design parameters are non-zero
int lsinc,nsinc,nsincm1,ishift;
double dsinc;
float[][] asinc;
}
/**
* Builds a table of interpolators for specified design parameters.
* Exactly one of the design parameters must be zero.
*/
private static Table makeTable(Design design) {
double emax = design.emax;
double fmax = design.fmax;
int lmax = design.lmax;
// The Kaiser window transition width is twice the difference
// between the Nyquist frequency 0.5 and the maximum frequency.
double wwin = 2.0*(0.5-fmax);
// The Kaiser window accounts for a hard-wired fraction of the maximum
// interpolation error. The other error will be due to table lookup.
double ewin = emax*EWIN_FRAC;
KaiserWindow kwin;
// If maximum frequency and length are specified, compute the error
// due to windowing. That windowing error is three times the error
// reported by the Kaiser window, because the conventional Kaiser
// window stopband error is aliased multiple times with the passband
// error. (A factor of at least two is necessary to handle the first
// aliasing, but experiments showed this factor to be inadequate. The
// factor three was found in testing to be sufficient for a wide range
// of design parameters.) Also, apply a lower bound to the error, so
// that the table of sinc approximations does not become too large.
if (emax==0.0) {
kwin = KaiserWindow.fromWidthAndLength(wwin,lmax);
ewin = 3.0*kwin.getError();
emax = ewin/EWIN_FRAC;
double etabMin = 1.1*PI*fmax/(NTAB_MAX-1);
double emaxMin = etabMin/(1.0-EWIN_FRAC);
if (emax0.0)?etab/(PI*fmax):1.0;
int nsincMin = 1+(int)ceil(1.0/dsinc);
int nsinc = 2;
while (nsinc _tables = new HashMap();
private static Table getTable(double emax, double fmax, int lmax) {
Design design = new Design(emax,fmax,lmax);
synchronized(_tables) {
Table table = _tables.get(design);
if (table==null)
table = makeTable(design);
return table;
}
}
private float interpolate(
double xscale, double xshift, int nxum, int nxu,
float[] yu, double x)
{
// Which uniform samples?
double xn = xshift+x*xscale;
int ixn = (int)xn;
int kyu = _ishift+ixn;
// Which sinc approximation?
double frac = xn-ixn;
if (frac<0.0)
frac += 1.0;
int ksinc = (int)(frac*_nsincm1+0.5);
float[] asinc = _asinc[ksinc];
// If no extrapolation is necessary, use a fast loop.
// Otherwise, extrapolate uniform samples, as necessary.
float yr = 0.0f;
if (kyu>=0 && kyu<=nxum) {
for (int isinc=0; isinc<_lsinc; ++isinc,++kyu)
yr += yu[kyu]*asinc[isinc];
} else if (_extrap==Extrapolation.ZERO) {
for (int isinc=0; isinc<_lsinc; ++isinc,++kyu) {
if (0<=kyu && kyu=0 && kyu<=nxum) {
for (int isinc=0; isinc<_lsinc; ++isinc,++kyu)
yu[kyu] += y*asinc[isinc];
} else if (_extrap==Extrapolation.ZERO) {
for (int isinc=0; isinc<_lsinc; ++isinc,++kyu) {
if (0<=kyu && kyu=0 && ky1u<=nx1um && ky2u>=0 && ky2u<=nx2um) {
for (int i2sinc=0; i2sinc<_lsinc; ++i2sinc,++ky2u) {
float asinc22 = asinc2[i2sinc];
float[] yuk2 = yu[ky2u];
float yr2 = 0.0f;
for (int i1sinc=0,my1u=ky1u; i1sinc<_lsinc; ++i1sinc,++my1u)
yr2 += yuk2[my1u]*asinc1[i1sinc];
yr += asinc22*yr2;
}
} else if (_extrap==Extrapolation.ZERO) {
for (int i2sinc=0; i2sinc<_lsinc; ++i2sinc,++ky2u) {
if (0<=ky2u && ky2u=0 && ky1u<=nx1um &&
ky2u>=0 && ky2u<=nx2um &&
ky3u>=0 && ky3u<=nx3um) {
for (int i3sinc=0; i3sinc<_lsinc; ++i3sinc,++ky3u) {
float asinc33 = asinc3[i3sinc];
float[][] yu3 = yu[ky3u];
float yr2 = 0.0f;
for (int i2sinc=0,my2u=ky2u; i2sinc<_lsinc; ++i2sinc,++my2u) {
float asinc22 = asinc2[i2sinc];
float[] yu32 = yu3[my2u];
float yr1 = 0.0f;
for (int i1sinc=0,my1u=ky1u; i1sinc<_lsinc; ++i1sinc,++my1u)
yr1 += yu32[my1u]*asinc1[i1sinc];
yr2 += asinc22*yr1;
}
yr += asinc33*yr2;
}
} else if (_extrap==Extrapolation.ZERO) {
for (int i3sinc=0; i3sinc<_lsinc; ++i3sinc,++ky3u) {
if (0<=ky3u && ky3u=0 && kyu<=nxum) {
for (int isinc=0; isinc<_lsinc; ++isinc,++kyu) {
int jyu = 2*kyu;
float asinci = asinc[isinc];
yr += yu[jyu ]*asinci;
yi += yu[jyu+1]*asinci;
}
} else if (_extrap==Extrapolation.ZERO) {
for (int isinc=0; isinc<_lsinc; ++isinc,++kyu) {
if (0<=kyu && kyu