Compute DCT using FFT
171776
def dct(data):
    
    N = len(data)//2
    fftdata     = fftpack.fft(data, axis=0)[:N+1]
    fftdata     /= N
    fftdata[0]  /= 2.
    fftdata[-1] /= 2.
    if np.isrealobj(data):
        data = np.real(fftdata)
    else:
        data = fftdata
    return data
46863
public static void FFT(ComplexNumber[] data, Direction direction) {
        double[] real = ComplexNumber.getReal(data);
        double[] img = ComplexNumber.getImaginary(data);
        if (direction == Direction.Forward)
            FFT(real, img);
        else
            FFT(img, real);
        if (direction == Direction.Forward) {
            for (int i = 0; i < real.length; i++) {
                data[i] = new ComplexNumber(real[i], img[i]);
            }
        } else {
            int n = real.length;
            for (int i = 0; i < n; i++) {
                data[i] = new ComplexNumber(real[i] / n, img[i] / n);
            }
        }
    }
923658
def _calc_fourier_spectrum(self, fa_length=None):
        """"""

        if fa_length is None:
            # Use the next power of 2 for the length
            n = 1
            while n < self.accels.size:
                n <<= 1
        else:
            n = fa_length

        self._fourier_amps = np.fft.rfft(self._accels, n)

        freq_step = 1. / (2 * self._time_step * (n / 2))
        self._freqs = freq_step * np.arange(1 + n / 2)
711559
def plot_fft(f, S, dt):
    '''
    '''
    import numpy

    xf = numpy.linspace(0.0, 1/(2.0*dt), N/2)
    plt.plot(xf, 2.0/N * numpy.abs(S[:N//2]), linewidth=_linewidth)
    plt.show()

    return None
39414
@Override
	public void compute(int x, int y, Point2D_F64 out ) {
		// grab precomputed normalized image coordinate at canonical location
		Point3D_F64 v = vectors[y*outWidth+x];
		// move to requested orientation
		GeometryMath_F64.mult(R,v,n); // TODO make faster by not using an array based matrix

		// compute pixel coordinate
		tools.normToEquiFV(n.x,n.y,n.z,out);
	}
862749
protected static void computeCosineWindow( GrayF64 cosine ) {
		double cosX[] = new double[ cosine.width ];
		for( int x = 0; x < cosine.width; x++ ) {
			cosX[x] = 0.5*(1 - Math.cos( 2.0*Math.PI*x/(cosine.width-1) ));
		}
		for( int y = 0; y < cosine.height; y++ ) {
			int index = cosine.startIndex + y*cosine.stride;
			double cosY = 0.5*(1 - Math.cos( 2.0*Math.PI*y/(cosine.height-1) ));
			for( int x = 0; x < cosine.width; x++ ) {
				cosine.data[index++] = cosX[x]*cosY;
			}
		}
	}
19555
def fft_plan(shape, dtype=np.complex64, axes=None, fast_math=True):
    
    # if not axes is None and any([a<0 for a in axes]):
    #     raise NotImplementedError("indices of axes have to be non negative, but are: %s"%str(axes))

    axes = _convert_axes_to_absolute(shape, axes)

    mock_buffer = MockBuffer(dtype, shape)

    fft_plan = FFT(mock_buffer, axes=axes).compile(cluda.ocl_api().Thread(get_device().queue),
                                                   fast_math=fast_math)

    return fft_plan
257987
def fft_freqs(n_fft, fs):
    
    return np.arange(0, (n_fft // 2 + 1)) / float(n_fft) * float(fs)
867605
public static void FFT2(ComplexNumber[][] data, Direction direction) {
        int n = data.length;
        int m = data[0].length;
        //ComplexNumber[] row = new ComplexNumber[m];//Math.max(m, n)];

        for (int i = 0; i < n; i++) {
            // copy row
            //for ( int j = 0; j < m; j++ )
            //row[j] = data[i][j];
            ComplexNumber[] row = data[i];
            // transform it
            FourierTransform.FFT(row, direction);
            // copy back
            for (int j = 0; j < m; j++)
                data[i][j] = row[j];
        }

        // process columns
        ComplexNumber[] col = new ComplexNumber[n];

        for (int j = 0; j < m; j++) {
            // copy column
            for (int i = 0; i < n; i++)
                col[i] = data[i][j];
            // transform it
            FourierTransform.FFT(col, direction);
            // copy back
            for (int i = 0; i < n; i++)
                data[i][j] = col[i];
        }
    }
40356
def fourier(x, N):
    """"""

    term = 0.
    for n in range(1, N, 2):
        term += (1. / n) * math.sin(n * math.pi * x / L)
    return (4. / (math.pi)) * term
237248
public void transform (double data[], int i0, int stride) {
    checkData(data,i0,stride);
    int p, p_1, q;

    if (n == 1) return;		/* identity operation */

    /* bit reverse the ordering of input data for decimation in time algorithm */
    bitreverse(data, i0, stride);

    /* apply fft recursion */
    p = 1; q = n ;
    for (int i = 1; i <= logn; i++) {
      int a, b;

      p_1 = p ;
      p = 2 * p ;
      q = q / 2 ;

      /* a = 0 */

      for (b = 0; b < q; b++) {
	double t0_real = data[i0+stride*b*p] + data[i0+stride*(b*p + p_1)];
	double t1_real = data[i0+stride*b*p] - data[i0+stride*(b*p + p_1)];
	  
	data[i0+stride*b*p] = t0_real;
	data[i0+stride*(b*p + p_1)] = t1_real;
      }

      /* a = 1 ... p_{i-1}/2 - 1 */

      {
	double w_real = 1.0;
	double w_imag = 0.0;

	double theta = - 2.0 * Math.PI / p;
	
	double s = Math.sin(theta);
	double t = Math.sin(theta / 2.0);
	double s2 = 2.0 * t * t;
	
	for (a = 1; a < (p_1)/2; a++) {
	  /* trignometric recurrence for w-> exp(i theta) w */
	    
	  {
	    double tmp_real = w_real - s * w_imag - s2 * w_real;
	    double tmp_imag = w_imag + s * w_real - s2 * w_imag;
	    w_real = tmp_real;
	    w_imag = tmp_imag;
	  }
	    
	  for (b = 0; b < q; b++) {
	    double z0_real = data[i0+stride*(b*p + a)];
	    double z0_imag = data[i0+stride*(b*p + p_1 - a)];
	    double z1_real = data[i0+stride*(b*p + p_1 + a)];
	    double z1_imag = data[i0+stride*(b*p + p - a)];
		
	    /* t0 = z0 + w * z1 */
	    data[i0+stride*(b*p + a)]    =z0_real + w_real * z1_real - w_imag * z1_imag;
	    data[i0+stride*(b*p + p - a)]=z0_imag + w_real * z1_imag + w_imag * z1_real;
	    /* t1 = -(z0 - w * z1) */
	    data[i0+stride*(b*p + p_1 - a)]=z0_real - w_real * z1_real + w_imag * z1_imag;
	    data[i0+stride*(b*p + p_1 + a)]=-(z0_imag - w_real * z1_imag - w_imag * z1_real);
	  }
	}
      }

      if (p_1 >  1) {
	for (b = 0; b < q; b++) {
	  /* a = p_{i-1}/2 */
	  data[i0+stride*(b*p + p - p_1/2)] *= -1 ;
	}}
    }
  }
272271
def dctmat(N,K,freqstep,orthogonalize=True):
    """"""
    cosmat = numpy.zeros((N, K), 'double')
    for n in range(0,N):
        for k in range(0, K):
            cosmat[n,k] = numpy.cos(freqstep * (n + 0.5) * k)
    if orthogonalize:
        cosmat[:,0] = cosmat[:,0] * 1./numpy.sqrt(2)
    return cosmat
1007038
public void computePixel(int which, double x, double y, Point2D_F64 output) {
		SurfaceRect r = scene.get(which);

		Point3D_F64 p3 = new Point3D_F64(-x,-y,0);
		SePointOps_F64.transform(r.rectToCamera, p3, p3);

		// unit sphere
		p3.scale(1.0/p3.norm());

		sphereToPixel.compute(p3.x,p3.y,p3.z,output);
	}
777449
def calc_spectrum(signal, rate):
    """"""
    npts = len(signal)
    padto = 1 << (npts - 1).bit_length()
    # print 'length of signal {}, pad to {}'.format(npts, padto)
    npts = padto

    sp = np.fft.rfft(signal, n=padto) / npts
    # print('sp len ', len(sp))
    freq = np.arange((npts / 2) + 1) / (npts / rate)
    # print('freq len ', len(freq))
    return freq, abs(sp)
1156345
private static void transformRadix2(double[] real, double[] imag) {
        int n = real.length;
        int levels = 31 - Integer.numberOfLeadingZeros(n);  // Equal to floor(log2(n))
//        if (1 << levels != n)
//            throw new IllegalArgumentException("Length is not a power of 2");
        double[] cosTable = new double[n / 2];
        double[] sinTable = new double[n / 2];
        for (int i = 0; i < n / 2; i++) {
            cosTable[i] = Math.cos(2 * Math.PI * i / n);
            sinTable[i] = Math.sin(2 * Math.PI * i / n);
        }

        // Bit-reversed addressing permutation
        for (int i = 0; i < n; i++) {
            int j = Integer.reverse(i) >>> (32 - levels);
            if (j > i) {
                double temp = real[i];
                real[i] = real[j];
                real[j] = temp;
                temp = imag[i];
                imag[i] = imag[j];
                imag[j] = temp;
            }
        }

        // Cooley-Tukey decimation-in-time radix-2 FFT
        for (int size = 2; size <= n; size *= 2) {
            int halfsize = size / 2;
            int tablestep = n / size;
            for (int i = 0; i < n; i += size) {
                for (int j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
                    double tpre = real[j + halfsize] * cosTable[k] + imag[j + halfsize] * sinTable[k];
                    double tpim = -real[j + halfsize] * sinTable[k] + imag[j + halfsize] * cosTable[k];
                    real[j + halfsize] = real[j] - tpre;
                    imag[j + halfsize] = imag[j] - tpim;
                    real[j] += tpre;
                    imag[j] += tpim;
                }
            }

            // Prevent overflow in 'size *= 2'
            if (size == n)
                break;
        }
    }
407639
def multiply_frequencies(signal, fs, frange, calibration_frequencies, attendB):
    """"""
    npts = len(signal)
    padto = 1 << (npts - 1).bit_length()

    X = np.fft.rfft(signal, n=padto)

    npts = padto
    f = np.arange((npts / 2) + 1) / (npts / fs)

    fidx_low = (np.abs(f - frange[0])).argmin()
    fidx_high = (np.abs(f - frange[1])).argmin()

    cal_func = interp1d(calibration_frequencies, attendB)
    roi = f[fidx_low:fidx_high]
    Hroi = cal_func(roi)
    H = np.zeros((len(X),))
    H[fidx_low:fidx_high] = Hroi

    H = smooth(H)
    # print 'H dB max', np.amax(H)

    H = 10 ** ((H).astype(float) / 20)
    # print 'H amp max', np.amax(H)

    # Xadjusted = X.copy()
    # Xadjusted[fidx_low:fidx_high] *= H
    # Xadjusted = smooth(Xadjusted)

    Xadjusted = X * H
    # print 'X max', np.amax(abs(X))
    # print 'Xadjusted max', np.amax(abs(Xadjusted))

    signal_calibrated = np.fft.irfft(Xadjusted)
    return signal_calibrated[:len(signal)]
592019
def fourier2dpad(data, zero_pad=True):
    
    if zero_pad:
        # zero padding size is next order of 2
        (N, M) = data.shape
        order = np.int(max(64., 2**np.ceil(np.log(2 * max(N, M)) / np.log(2))))

        # this is faster than np.pad
        datapad = np.zeros((order, order), dtype=float)
        datapad[:data.shape[0], :data.shape[1]] = data
    else:
        datapad = data

    # Fourier transform
    fft = np.fft.fftshift(np.fft.fft2(datapad))
    return fft