All Downloads are FREE. Search and download functionalities are using the official Maven repository.

com.aparapi.examples.extension.FFTExample Maven / Gradle / Ivy

/**
 * This product currently only contains code developed by authors
 * of specific components, as identified by the source code files.
 *
 * Since product implements StAX API, it has dependencies to StAX API
 * classes.
 *
 * For additional credits (generally to people who reported problems)
 * see CREDITS file.
 */
package com.aparapi.examples.extension;

import com.aparapi.*;
import com.aparapi.device.*;
import com.aparapi.internal.kernel.*;
import com.aparapi.opencl.*;
import com.aparapi.opencl.OpenCL.*;

import java.util.*;

public class FFTExample{

   @Resource("com/aparapi/examples/extension/fft.cl") interface FFT extends OpenCL{

      public FFT forward(//
            Range _range,//
            @GlobalReadWrite("real") float[] real,//
            @GlobalReadWrite("imaginary") float[] imaginary//
      );
   }

   static void fft(float[] x, float[] y) {
      final short dir = 1;
      final long m = 10;
      int n, i, i1, j, k, i2, l, l1, l2;
      double c1, c2, tx, ty, t1, t2, u1, u2, z;

      /* Calculate the number of points */
      n = 1;
      for (i = 0; i < m; i++) {
         n *= 2;
      }

      /* Do the bit reversal */
      i2 = n >> 1;
      j = 0;
      for (i = 0; i < (n - 1); i++) {
         if (i < j) {
            tx = x[i];
            ty = y[i];
            x[i] = x[j];
            y[i] = y[j];
            x[j] = (float) tx;
            y[j] = (float) ty;
         }
         k = i2;
         while (k <= j) {
            j -= k;
            k >>= 1;
         }
         j += k;
      }

      /* Compute the FFT */
      c1 = -1.0;
      c2 = 0.0;
      l2 = 1;
      for (l = 0; l < m; l++) {
         l1 = l2;
         l2 <<= 1;
         u1 = 1.0;
         u2 = 0.0;
         for (j = 0; j < l1; j++) {
            for (i = j; i < n; i += l2) {
               i1 = i + l1;
               t1 = (u1 * x[i1]) - (u2 * y[i1]);
               t2 = (u1 * y[i1]) + (u2 * x[i1]);
               x[i1] = (float) (x[i] - t1);
               y[i1] = (float) (y[i] - t2);
               x[i] += (float) t1;
               y[i] += (float) t2;
            }
            z = (u1 * c1) - (u2 * c2);
            u2 = (u1 * c2) + (u2 * c1);
            u1 = z;
         }
         c2 = Math.sqrt((1.0 - c1) / 2.0);
         if (dir == 1) {
            c2 = -c2;
         }
         c1 = Math.sqrt((1.0 + c1) / 2.0);
      }

      /* Scaling for forward transform */
      /*if (dir == 1) {
         for (i=0;i 0.01) {
            System.out.printf("%d %5.2f %5.2f %5.2f\n", i, initial[i], real[i], referenceReal[i]);
         }
      }

   }
}




© 2015 - 2024 Weber Informatics LLC | Privacy Policy