/************Chrip Z Transform *****************
   * N = # input samples. M = # output samples.
   *
   * fsam   = the sample frequency in Hz.
   * fstart = the start frequency in Hz.
   * fstop  = the stop frequency in Hz for the band over which the transform is computed.
   *
   * (fstart - fstop)/M = new resolution  
   *
   * Note: this method returns an array of length L. 
   * L = the returned transform length. L will always be larger than M.
   * See code for how L is determined   
   *
   **********************************************/
  public double[ ][ ] czt( double [ ][ ] array, int N, int M, double fStart, double fStop, double fSam )
  {
    int L;

    /***Determine length of CZT output***/
    if( (N+M) > 512){
      L = 1024;
    }else if( ((N+M) > 256) && ((N+M) <= 512) ){
      L = 512;
    }else if( ((N+M) > 128) && ((N+M) <= 256) ){
      L = 256;
    }else if( ((N+M) > 64 ) && ((N+M) <= 128) ){
      L = 128;
    }else if( ((N+M) > 32 ) && ((N+M) <= 64) ){
      L = 64;
    }else if( ((N+M) > 16 ) && ((N+M) <= 32) ){
      L = 32;
    }else{
      L = 16;
    }

    double[ ][ ] g = new double[L][2];
    double[ ][ ] h = new double[L][2];
    double theta0;
    double phi0;
    double psi;
    double a;
    double b;
    int n;
    int k;

    phi0 = 2.*Math.PI*(fStop-fStart)/fSam/(M-1);
    theta0 = 2.*Math.PI*fStart/fSam;

    /*** Create arc coefficients ***/
    for( n = 0; n < M; n++ ){
      h[n][0] = Math.cos( n*n/2.*phi0 );
      h[n][1] = Math.sin( n*n/2.*phi0 );
    }
    for( n = M; n < L-N; n++ ){
      h[n][0] = 0.;
      h[n][1] = 0.;
    }
    for( n = L-N; n < L; n++){
      h[n][0] = Math.cos( (L-n)*(L-n)/2.*phi0 );
      h[n][1] = Math.sin( (L-n)*(L-n)/2.*phi0 );
    }

    /*** Prepare signal ***/
    for( n = 0; n < N; n++ ){
      g[n][0] = array[n][0];
      g[n][1] = array[n][1];
    }
    for( n = N; n < L; n++ ){
      g[n][0] = 0.;
      g[n][1] = 0.;
    }
    
    double s;
    double c;
    
    for( n = 0; n < N; n++ ){
      psi = n*theta0 + n*n/2.*phi0;
      c =  Math.cos( psi);
      s = -Math.sin( psi );
      a = c*g[n][0] - s*g[n][1];
      b = s*g[n][0] + c*g[n][1];
      g[n][0] = a;
      g[n][1] = b;
    }
                       //use your favorite fft algorithm here.
    g = fft_1d( g );   //fft of samples    
    h = fft_1d( h );   //fft of arc coeff
    
    /** convolution in the time domain is multiplication in the frequency domain **/
    /** multiplication in the time domain is convolution in the frequency domain **/
    for( n = 0; n < L; n++ ){
      c = g[n][0];
      s = g[n][1];
      a = c*h[n][0] - s*h[n][1];
      b = s*h[n][0] + c*h[n][1];
      g[n][0] = a/L;  /* for scaling purposes since
			 fft_1d does not use scale */
      g[n][1] = b/L;
    }

    g = ifft_1d( g );   //use your favorite fft algorithm here.

    for( k = 0; k < M; k++ ){
      psi = k*k/2.*phi0;
      c =  Math.cos( psi );
      s = -Math.sin( psi );
      a = c*g[k][0] - s*g[k][1];
      b = s*g[k][0] + c*g[k][1];
      g[k][0] = a;
      g[k][1] = b;
    }
    
    return g;
  
  }