org.apache.commons.math3.transform.TransformType Java Examples

/**
  * Forward or Inverse *squared* fast fourier transform onto a particular axis
  *
* N.B. returns the squared complex fourier component!
*
  * @param axis      the axis (0 is rows, 1 is columns)
  * @param direction TransformType.FORWARD or TransformType.INVERSE
  * @return The new matrix with a fft or ifft applied to each row or column
  */
 public Matrix fft2(int axis, TransformType direction) {
     ParameterChecker.checkAxis(axis);

     FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
     double[][] ft = new double[this.getRowDimension()][this.getColumnDimension()];
     if (axis == 0) {
         for (int c = 0; c < this.getColumnDimension(); c++) {
             Complex[] complexResult = fft.transform(this.getColumn(c), direction);
             for (int i = 0; i < complexResult.length; i++)
                 ft[i][c] = complexResult[i].abs()*complexResult[i].abs();
         }
     } else {
	// TODO: This is inefficient....
         return this.transpose().fft2(0, direction).transpose();
     }
     return new Matrix(ft);
 }
 

@Override
public Object doWork(Object v) throws IOException {

  double[] data = ((List<?>)v).stream().mapToDouble(value -> ((Number)value).doubleValue()).toArray();

  FastFourierTransformer fastFourierTransformer = new FastFourierTransformer(DftNormalization.STANDARD);
  Complex[] complex = fastFourierTransformer.transform(data, TransformType.FORWARD);

  double[] real = new double[complex.length];
  double[] imaginary = new double[complex.length];

  for(int i=0; i<real.length; ++i) {
    real[i] = complex[i].getReal();
    imaginary[i] = complex[i].getImaginary();
  }

  double[][] d = new double[2][];
  d[0]=real;
  d[1]=imaginary;

  Matrix matrix = new Matrix(d);
  matrix.setRowLabels(clabels);
  return matrix;
}
 

@Override
public void consume(List<Datum> records) throws Exception {
    for (Datum d: records){
        metricVector = d.metrics();
        // TODO: look for decent FFT implementation that doesn't need pwr of 2
        nextPowTwo = Math.max(2,2*Integer.highestOneBit(metricVector.getDimension()-1));
        paddedInput = metricVector.append(new ArrayRealVector(nextPowTwo - metricVector.getDimension()));

        transformer = new FastFourierTransformer(DftNormalization.STANDARD);
        FFTOutput = transformer.transform(paddedInput.toArray(), TransformType.FORWARD);
        transformedMetricVector = new ArrayRealVector();
        for (Complex c: FFTOutput){
            transformedMetricVector = transformedMetricVector.append(c.getReal());
            transformedMetricVector = transformedMetricVector.append(c.getImaginary());
        }
        output.add(new Datum(d, transformedMetricVector));
    }
}
 

/**
 *
 * @param shot
 */
@Override
public void processSegment(SegmentContainer shot) {

    BufferedImage image = shot.getAvgImg().getBufferedImage();

    List<Contour> contours = ContourHelper.getContours(image);
    List<Point2D_I32> contour = contours.get(0).internal.get(0);

    if (image != null) {
        FastFourierTransformer transformer = new FastFourierTransformer(DftNormalization.STANDARD);
        double[] distancefunction = ContourHelper.centroidDistance(contour, true);
        Complex[] signature = transformer.transform(distancefunction, TransformType.FORWARD);
        float[] descriptors = new float[DESCRIPTOR_LENGTH];
        for (int i = 1;i<DESCRIPTOR_LENGTH;i++) {
            descriptors[i] = (float) (signature[i].abs() / signature[0].abs());
        }
        this.persist(shot.getId(), new FloatVectorImpl(descriptors));
    }
}
 

/**
 *
 * @param sc
 * @param qc
 * @return
 */
@Override
public List<ScoreElement> getSimilar(SegmentContainer sc, ReadableQueryConfig qc) {

    BufferedImage image = sc.getAvgImg().getBufferedImage();

    qc = setQueryConfig(qc);

    List<Contour> contours = ContourHelper.getContours(image);
    List<Point2D_I32> contour =  contours.get(0).internal.get(0);

    if (image != null) {
        FastFourierTransformer transformer = new FastFourierTransformer(DftNormalization.STANDARD);
        double[] distancefunction = ContourHelper.centroidDistance(contour, true);
        Complex[] signature = transformer.transform(distancefunction, TransformType.FORWARD);
        float[] descriptors = new float[DESCRIPTOR_LENGTH];
        for (int i = 1;i<DESCRIPTOR_LENGTH;i++) {
            descriptors[i] = (float) (signature[i].abs() / signature[0].abs());
        }
        return this.getSimilar(descriptors, qc);
    } else {
        return new ArrayList<>();
    }
}
 

public void evaluate(List<Datum> data) {
    double[] values = formatData(data);
    // FFT
    Complex[] fft = fftTran.transform(values, TransformType.FORWARD);
    // Multiply by complex conjugate
    for (int i = 0; i < fft.length; i ++) {
        fft[i] = fft[i].multiply(fft[i].conjugate());
    }
    // Inverse transform
    fft = fftTran.transform(fft, TransformType.INVERSE);

    correlations = new double[maxLag];
    for (int i = 1; i < maxLag; i++) {
        correlations[i] = fft[i].getReal() / fft[0].getReal();
    }
}
 

/**
 * Performs a forward fourier transformation on the provided, real valued data. The method makes sure,
 * that the size of the array is a power of two (for which the FFT class has been optimized) and pads
 * the data with zeros if necessary. Furthermore, one can provide a WindowingFunction that will be applied
 * on the data.
 *
 * <strong>Important: </strong>Every call to forward() replaces all the existing data in the current instance. I.e.
 * the same instance of FFT can be re-used.
 *
 * @param data Data to be transformed.
 * @param window WindowFunction to use for the transformation.
 */
public void forward(double[] data, float samplingrate, WindowFunction window) {
    this.windowFunction = window;
    this.samplingrate = samplingrate;

    int actual_length = data.length;
    int valid_length = FFTUtil.nextPowerOf2(actual_length);
    double[] reshaped = new double[valid_length];
    for (int i = 0; i<reshaped.length; i++) {
        if (i < actual_length) {
            reshaped[i] = data[i] * this.windowFunction.value(i, valid_length);
        } else {
            reshaped[i] = 0;
        }
    }

    /* Perform FFT using FastFourierTransformer library. */
    FastFourierTransformer transformer = new FastFourierTransformer(DftNormalization.STANDARD);
    this.data = transformer.transform(reshaped, TransformType.FORWARD);

    /* Reset the calculated properties. */
    this.powerSpectrum = null;
    this.magnitudeSpectrum = null;
}
 

public static double[] fft(short[] buffer, int offset, int len) {
        int len2 = (int) Math.pow(2, Math.ceil(Math.log(len) / Math.log(2)));

        final double[][] dataRI = new double[][]{
                new double[len2], new double[len2]
        };

        double[] dataR = dataRI[0];
        double[] dataI = dataRI[1];

        double powerInput = 0;
        for (int i = 0; i < len; i++) {
            dataR[i] = buffer[offset + i] / (float) 0x7fff;
            powerInput += dataR[i] * dataR[i];
        }
        powerInput = Math.sqrt(powerInput / len);

        FastFourierTransformer.transformInPlace(dataRI, DftNormalization.STANDARD, TransformType.FORWARD);

        double[] data = new double[len2 / 2];

        data[0] = 10 * Math.log10(Math.pow(new Complex(dataR[0], dataI[0]).abs() / len2, 2));

        double powerOutput = 0;
        for (int i = 1; i < data.length; i++) {
            Complex c = new Complex(dataR[i], dataI[i]);
            double p = c.abs();
            p = p / len2;
            p = p * p;
            p = p * 2;
            double dB = 10 * Math.log10(p);

            powerOutput += p;
            data[i] = dB;
        }
        powerOutput = Math.sqrt(powerOutput);

//        if(powerInput != powerOutput) {
//            throw new RuntimeException("in " + powerInput + " out " + powerOutput);
//        }

        return data;
    }
 

/**
 * Method to generate a frequency spectrum of the frame using FFT
 * 
 * @param frame
 *            Frame to be transformed
 * @param Fs
 *            Sampling Frequency
 * @return an Array showing the realtive powers of all frequencies
 */
private static double[] transformFrameFFT(double[] frame, int Fs) {
	final FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
	final Complex[] spectrum = fft.transform(frame, TransformType.FORWARD);
	final double[] powerSpectrum = new double[frame.length / 2 + 1];
	for (int ii = 0; ii < powerSpectrum.length; ii++) {
		final double abs = spectrum[ii].abs();
		powerSpectrum[ii] = abs * abs;
	}
	return powerSpectrum;
}
 

/**
 * Extracts the Lightfield Fourier descriptors from a provided BufferedImage. The returned list contains
 * elements for each identified contour of adequate size.
 *
 * @param image Image for which to extract the Lightfield Fourier descriptors.
 * @param poseidx Poseidx of the extracted image.
 * @return List of descriptors for image.
 */
@Override
protected List<float[]> featureVectorsFromImage(BufferedImage image, int poseidx) {
    final List<Contour> contours = ContourHelper.getContours(image);
    final List<float[]> features = new ArrayList<>();

    /* Select the largest, inner contour from the list of available contours. */
    for (Contour contour : contours) {
        for (List<Point2D_I32> inner : contour.internal) {
            /* Check size of selected contour. */
            if (inner.size() < SIZE * 2) {
              continue;
            }

            /* Calculate the descriptor for the selected contour. */
            double[] cds = ContourHelper.centroidDistance(inner, true);
            Complex[] results = this.transformer.transform(cds, TransformType.FORWARD);
            double magnitude = results[0].abs();
            float[] feature = new float[SIZE];
            for (int i = 1; i < SIZE; i++) {
                feature[i] = (float) (results[i+1].abs() / magnitude);
            }
            feature = MathHelper.normalizeL2InPlace(feature);
            feature[0] = poseidx;
            features.add(feature);
        }
    }

    return features;
}
 

@Override
public Object doWork(Object v) throws IOException {

  if(v instanceof Matrix) {

    Matrix matrix = (Matrix)v;
    double[][] data = matrix.getData();
    double[] real = data[0];
    double[] imaginary = data[1];
    Complex[] complex = new Complex[real.length];

    for (int i = 0; i < real.length; ++i) {
     complex[i] = new Complex(real[i], imaginary[i]);
    }

    FastFourierTransformer fastFourierTransformer = new FastFourierTransformer(DftNormalization.STANDARD);
    Complex[] result  = fastFourierTransformer.transform(complex, TransformType.INVERSE);

    List<Number> realResult = new ArrayList<>();
    for (int i = 0; i < result.length; ++i) {
      realResult.add(result[i].getReal());
    }

    return realResult;
  } else {
    throw new IOException("ifft function requires a matrix as a parameter");
  }
}
 

public static void main(String[] args) {
	final double Fs = 8000;
	final int N = 256;

	// *** Generate a signal with frequency f
	double[] signal = new double[N];
	System.out.println("Frame length: " + N + " samples => " + (N / Fs) + "ms");

	final double f = 1234.0; // Hz
	for (int ii = 0; ii < N; ii++) {
		signal[ii] = Math.sin(2.0 * Math.PI * f * (double) ii / Fs);
	}

	// Should apply hamming window to the frame.

	// *** Get power spectrum
	final FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);

	// Note: signal.length should have been a power of 2
	final Complex[] spectrum = fft.transform(signal, TransformType.FORWARD);
	final double[] powerSpectrum = new double[N / 2 + 1];
	for (int ii = 0; ii < powerSpectrum.length; ii++) {
		final double abs = spectrum[ii].abs();
		powerSpectrum[ii] = abs * abs;
	}

	// Expect a sharp peak at around frequency f - index f/Fs *
	// signal.length
	final double center = (f / Fs) * N;
	final int start = (int) (center - 10.0);
	final int end = (int) (center + 10.0);

	System.out.println("Center frequency in the FFT: " + center);
	for (int ii = start; ii < end; ii++) {
		System.out.format("% 3d (% 4.0fHz) power:%01.01f\n", ii, ii * Fs / N, powerSpectrum[ii]);
	}

	// *** Detect
	double totalPower = sum(powerSpectrum);
	double signalPower = powerSpectrum[39] + powerSpectrum[40];
	double ratio = signalPower / totalPower;

	// around 80%. Reason being leakage in the FFT. The actual spectrum is a
	// sigmoid function with sidelobes. If you use e.g. a Hamming window,
	// the side-lobes are suppressed at the cost an even broader center
	// lobe. It's ok though, it just affects the width of the band where you
	// look for the actual frequency of interest.
	System.out.println("Detection ratio: " + ratio);
}
 

public static void highpassfftfilt(List<double[]> set, double lowfreq, double highfreq) {
	/*
	 * steps:
	 * 1. pad data with 0s to get the number of samples to a power of 2
	 * 2. perform fft
	 * 3. filter based on frequency (see http://stackoverflow.com/a/2876292)
	 * 4. inverse fft
	 */

	//step 1
	ArrayList<double[]> padded = new ArrayList<double[]>(set);
	int padTo = Filters.findNextLargestPower2(padded.size());
	for(int i = padded.size(); i < padTo; i++) {
		padded.add(new double[]{0.0, 0.0});
	}
	double[][] arr = Filters.toArray(padded);

	//step 2
	FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
	Complex[] out = fft.transform(arr[1], TransformType.FORWARD);

	//step 3
	double sampFreq = 1000.0 / (arr[0][1]-arr[0][0]);
	int posHalf = out.length / 2;
	for(int i = 0; i < posHalf; i++) {
		int negInd = out.length - 1 - i;
		double currFreq = (double)i * (sampFreq / (double)posHalf);
		if (currFreq > lowfreq) {
			out[i] = Complex.ZERO;
			out[negInd] = Complex.ZERO;
		} /* else if (currFreq < highfreq) {
			double scale = 1.0 - ((currFreq - highFreq) / (lowFreq - highFreq));
			out[i] = out[i].multiply(scale);
			out[negInd] = out[negInd].multiply(scale);
		} */
	}

	//step 4
	out = fft.transform(out, TransformType.INVERSE);

	//write changes
	for(int i = 0; i < set.size(); i++) {
		set.set(i, new double[]{set.get(i)[0], out[i].getReal()});
	}
}
 

public Instances process(Instances input) throws Exception {
    Instances instances = new Instances(input);
    double[][] data = null;
    double[][] FFTData = null;
    Instances FFTInstances = null;
    FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
    Complex[] complexData = null;

    data = fromWekaInstancesArray(instances, true);
    FFTInstances = determineOutputFormat(instances);

    FFTData = new double[instances.size()][nfft/2 + 1];

    for (int i = 0; i < FFTData.length; i++){

        complexData = new Complex[nfft];
        for (int j = 0; j < complexData.length; j++){
            complexData[j] = new Complex(0.0, 0.0);
        }

        double mean = 0;
        if (data[i].length < nfft){
            for (int j = 0; j < data[i].length; j++) {
                mean += data[i][j];
            }
            mean /= data[i].length;
        }

        //int limit = nfft < data[i].length ? nfft : data[i].length;
        for (int j = 0; j < nfft; j++) {
            if(j < data[i].length)
                complexData[j] = new Complex(data[i][j], 0);
            else
                complexData[j] = new Complex(mean, 0);
        }

        complexData = fft.transform(complexData, TransformType.FORWARD);

        for (int j = 0; j < (nfft/2); j++) {
            FFTData[i][j] = complexData[j].abs();
        }

        FFTData[i][FFTData[i].length - 1] = instances.get(i).classValue();
    }

    for (int i = 0; i < FFTData.length; i++) {
        FFTInstances.add(new DenseInstance(1.0, FFTData[i]));
    }

    return FFTInstances;
}
 

@Override
public Map<String, Function<Double, Double>> getValue(final double evaluationTime, final CharacteristicFunctionModel model) throws CalculationException {

	final CharacteristicFunction modelCF = model.apply(getMaturity());

	final double lineOfIntegration = 0.5 * (getIntegrationDomainImagUpperBound()+getIntegrationDomainImagLowerBound());

	final double lambda = 2*Math.PI/(numberOfPoints*gridSpacing); //Equation 23 Carr and Madan
	final double upperBound = (numberOfPoints * lambda)/2.0; //Equation 20 Carr and Madan

	final Complex[] integrandEvaluations = new Complex[numberOfPoints];

	for(int i = 0; i<numberOfPoints; i++) {

		final double u = gridSpacing * i;

		//Integration over a line parallel to the real axis
		final Complex z = new Complex(u,-lineOfIntegration);

		//The characteristic function is already discounted
		final Complex numerator = modelCF.apply(z.subtract(Complex.I));

		final Complex denominator = apply(z);
		Complex ratio = numerator.divide(denominator);
		ratio = (ratio.multiply(((Complex.I).multiply(upperBound*u)).exp())).multiply(gridSpacing);

		double delta;
		if (i==0){
			delta=1.0;
		}else{
			delta = 0.0;
		}
		final double simpsonWeight = (3+Math.pow(-1,i+1)-delta)/3;

		integrandEvaluations[i] = ratio.multiply(simpsonWeight);
	}

	//Compute the FFT
	Complex[] transformedVector = new Complex[numberOfPoints];
	final FastFourierTransformer fft=new FastFourierTransformer(DftNormalization.STANDARD);
	transformedVector=	fft.transform(integrandEvaluations,TransformType.FORWARD);

	//Find relevant prices via interpolation
	final double[] logStrikeVector = new double[numberOfPoints];
	final double[] strikeVector = new double[numberOfPoints];
	final double[] optionPriceVector = new double[numberOfPoints];

	for(int j = 0; j<numberOfPoints; j++) {
		logStrikeVector[j] = -upperBound+lambda*j;
		strikeVector[j] = Math.exp(logStrikeVector[j]);
		optionPriceVector[j] = (transformedVector[j].multiply(Math.exp(-lineOfIntegration * logStrikeVector[j]))).getReal()/Math.PI;
	}

	final RationalFunctionInterpolation interpolation = new RationalFunctionInterpolation(strikeVector, optionPriceVector,intMethod, extMethod);

	final Complex minusI = new Complex(0,-1);
	final double residueTerm = (modelCF.apply(minusI)).getReal();

	final Function<Double, Double> strikeToPrice = new Function<Double, Double>(){

		@Override
		public Double apply(final Double t) {
			return residueTerm + interpolation.getValue(t);
		}

	};

	final HashMap<String, Function<Double, Double>> results = new HashMap<>();
	results.put("valuePerStrike", strikeToPrice);
	return results;
}
 

@Override
public Map<String, Function<Double, Double>> getValue(final double evaluationTime, final CharacteristicFunctionModel model) throws CalculationException {

	final CharacteristicFunction modelCF = model.apply(getMaturity());

	final double lineOfIntegration = 0.5 * (getIntegrationDomainImagUpperBound()+getIntegrationDomainImagLowerBound());

	final double lambda = 2*Math.PI/(numberOfPoints*gridSpacing); //Equation 23 Carr and Madan
	final double upperBound = (numberOfPoints * lambda)/2.0; //Equation 20 Carr and Madan

	final Complex[] integrandEvaluations = new Complex[numberOfPoints];

	for(int i = 0; i<numberOfPoints; i++) {

		final double u = gridSpacing * i;

		//Integration over a line parallel to the real axis
		final Complex z = new Complex(u,-lineOfIntegration);

		//The characteristic function is already discounted
		final Complex numerator = modelCF.apply(z.subtract(Complex.I));

		final Complex denominator = apply(z);
		Complex ratio = numerator.divide(denominator);
		ratio = (ratio.multiply(((Complex.I).multiply(upperBound*u)).exp())).multiply(gridSpacing);

		double delta;
		if (i==0){
			delta=1.0;
		}else{
			delta = 0.0;
		}
		final double simpsonWeight = (3+Math.pow(-1,i+1)-delta)/3;

		integrandEvaluations[i] = ratio.multiply(simpsonWeight);
	}

	//Compute the FFT
	Complex[] transformedVector = new Complex[numberOfPoints];
	final FastFourierTransformer fft=new FastFourierTransformer(DftNormalization.STANDARD);
	transformedVector=	fft.transform(integrandEvaluations,TransformType.FORWARD);

	//Find relevant prices via interpolation
	final double[] logStrikeVector = new double[numberOfPoints];
	final double[] strikeVector = new double[numberOfPoints];
	final double[] optionPriceVector = new double[numberOfPoints];

	for(int j = 0; j<numberOfPoints; j++) {
		logStrikeVector[j] = -upperBound+lambda*j;
		strikeVector[j] = Math.exp(logStrikeVector[j]);
		optionPriceVector[j] = (transformedVector[j].multiply(Math.exp(-lineOfIntegration * logStrikeVector[j]))).getReal()/Math.PI;
	}

	final RationalFunctionInterpolation interpolation = new RationalFunctionInterpolation(strikeVector, optionPriceVector,intMethod, extMethod);

	final Complex minusI = new Complex(0,-1);
	final double residueTerm = (modelCF.apply(minusI)).getReal();

	final Function<Double, Double> strikeToPrice = new Function<Double, Double>(){

		@Override
		public Double apply(final Double t) {
			return residueTerm + interpolation.getValue(t);
		}

	};

	final HashMap<String, Function<Double, Double>> results = new HashMap<>();
	results.put("valuePerStrike", strikeToPrice);
	return results;
}
 

/**
 * Forward fast fourier transform onto a particular axis
 *
 * @param axis the axis (0 is rows, 1 is columns)
 * @return The new matrix with a fft applied to each row or column
 */
public Matrix fft2(int axis) {
    return fft2(axis, TransformType.FORWARD);
}
 

/**
 * Inverse fast fourier transform onto a particular axis
 *
 * @param axis the axis (0 is rows, 1 is columns)
 * @return The new matrix with an ifft applied to each row or column
 */
public Matrix ifft2(int axis) {
    return fft2(axis, TransformType.INVERSE);
}