Package com.nr.sp

Source Code of com.nr.sp.DftInt

package com.nr.sp;

import org.netlib.util.doubleW;

import com.nr.fft.*;
import static java.lang.Math.*;
import com.nr.UniVarRealValueFun;
import com.nr.interp.Poly_interp;

/**
* Example program illustrating how to use the routine dftcor.
* The user supplies an external function func that returns the
* quantity h(t). The routine then returns S(a,b, cos(wt)*h(t)*dt as
* cosint and S(a,b, sin(wt)*h(t)*dt as sinint.
*
*/
public class DftInt {
  final int M=64,NDFT=1024,MPOL=6;
  /*
   The values of M, NDFT, and MPOL are merely illustrative and should be optimized for your
   particular application. M is the number of subintervals, NDFT is the length of the FFT ( a power
   of 2), and MPOL is the degree of polynomial interpolation used to obtain the desired frequency from the FFT
   */
  double aold = -1.e30,bold = -1.e30,delta;
  int init=0;
  double[] data = new double[NDFT];
  double[] endpts = new double[8];
  UniVarRealValueFun funcold;
 

  public void dftint(final UniVarRealValueFun func, final double a, final double b, final double w,
      final doubleW cosint, final doubleW sinint) {
    final double TWOPI=6.283185307179586476;
    int j,nn;
    double c,cdft,en,s,sdft;
    doubleW corfac = new doubleW(0);
    doubleW corim = new doubleW(0);
    doubleW corre = new doubleW(0);
    double[] cpol = new double[MPOL];
    double[] spol = new double[MPOL];
    double[] xpol = new double[MPOL];
    if (init != 1 || a != aold || b != bold || func != funcold) {
        // Do we need to initialize, or is only ω changed?
      init=1;
      aold=a;
      bold=b;
      funcold=func;
      delta=(b-a)/M;
      for (j=0;j<M+1;j++)   //    Load the function values into the data array.
        data[j]=func.funk(a+j*delta);
      for (j=M+1;j<NDFT;j++)   // Zero-pad the rest of the data array
        data[j]=0.0;
      for (j=0;j<4;j++) {    // Load the endpoints.
        endpts[j]=data[j];
        endpts[j+4]=data[M-3+j];
      }
      FFT.realft(data,1);  
      // realft returns the unused value corresponding to ω_{N/2} in data[1]. We actually want
      // this element to contain the imaginary part corresponding to ω_0, which is zero
      data[1]=0.0;
    }
    // Now interpolate on the DFT result for the desired frequency. If the frequency is an ω_n,
    // i.e. the quantity en is an integer, then cdft = data[2*en-2], sdft = ata[2*en-1],
    // and you could omit the interpolation
    en=w*delta*NDFT/TWOPI+1.0;
    nn=min(max((int)(en-0.5*MPOL+1.0),1),NDFT/2-MPOL+1)// Leftmost point for the interpolation
    for (j=0;j<MPOL;j++,nn++) {
      cpol[j]=data[2*nn-2];
      spol[j]=data[2*nn-1];
      xpol[j]=nn;
    }
    cdft = new Poly_interp(xpol,cpol,MPOL).interp(en);
    sdft = new Poly_interp(xpol,spol,MPOL).interp(en);
    Fourier.dftcor(w,delta,a,b,endpts,corre,corim,corfac);    // Now get the endpoint correction and the multiplicative
    cdft *= corfac.val;                                                         // factor W(θ)
    sdft *= corfac.val;
    cdft += corre.val;
    sdft += corim.val;
    c=delta*cos(w*a);       // Finally multiply by Δ and exp(ιωα)
    s=delta*sin(w*a);
    cosint.val=c*cdft-s*sdft;
    sinint.val=s*cdft+c*sdft;
  }
}
TOP

Related Classes of com.nr.sp.DftInt

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.