// Path: lib/src/main/java/hrv/lib/hrv/RRData.java
// public class RRData {
// 
// 	private final TimeUnit rrValueAxisUnit;
// 	private final /*@ spec_public */ TimeUnit rrTimeAxisUnit;
// 
// 	private final double[] rrValuesAxis;
// 	private final double[] rrTimeAxis;
// 
// 	public RRData(double[] timeAxis, TimeUnit timeAxisUnit, double[] valueAxis, TimeUnit valueAxisUnit) {
// 		this.rrValuesAxis = valueAxis;
// 		this.rrValueAxisUnit = valueAxisUnit;
// 
// 		this.rrTimeAxis = timeAxis;
// 		this.rrTimeAxisUnit = timeAxisUnit;
// 	}
// 
// 	/**
// 	 * Creates a new RRData object from an RR-Interval
// 	 * 
// 	 * @param rrInterval
// 	 *            RR-Interval
// 	 * @param unit
// 	 *            of the given RR-Data.
// 	 * @return new RRData object with given RR-Interval data.
// 	 */
// 	public static RRData createFromRRInterval(double[] rrInterval, TimeUnit unit) {
// 		var rrTimeAxis = new double[rrInterval.length];
// 		for (var i = 1; i < rrInterval.length; i++) {
// 			rrTimeAxis[i] = rrTimeAxis[i - 1] + rrInterval[i - 1];
// 		}
// 
// 		return new RRData(rrTimeAxis, unit, rrInterval, unit);
// 	}
// 
// 	public double[] getValueAxis() {
// 		return rrValuesAxis;
// 	}
// 
// 	public double[] getTimeAxis() {
// 		return rrTimeAxis;
// 	}
// 	
// 	public TimeUnit getValueAxisUnit() {
// 		return rrValueAxisUnit;
// 	}
// 
// 	public TimeUnit getTimeAxisUnit() {
// 		return rrTimeAxisUnit;
// 	}
// }
// 
// Path: lib/src/main/java/hrv/lib/hrv/calc/parameter/AvgSampleSizeCalculator.java
// public class AvgSampleSizeCalculator implements HRVDataProcessor {
// 
// 	@Override
// 	public HRVParameter process(RRData data) {
// 		double maxX = data.getTimeAxis()[data.getTimeAxis().length - 1];
// 		double minX = data.getTimeAxis()[0];
// 		double xLength = maxX - minX;
// 		return new HRVParameter(HRVParameterEnum.AVG_SAMPLE_SIZE, xLength / (data.getTimeAxis().length - 1),
// 				"1 / " + data.getTimeAxisUnit().toString());
// 	}
// 	
// 
// 	@Override
// 	public boolean validData(RRData data) {
// 		return data.getValueAxis().length > 0;
// 	}
// }


import hrv.lib.hrv.RRData;
import hrv.lib.hrv.calc.parameter.AvgSampleSizeCalculator;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.TransformType;
package hrv.lib.hrv.calc.psd;


/**
 * Uses FFT to estimate the power spectra density of RR-Data-Intervals
 * 
 * @author Julian
 *
 */
public class StandardPowerSpectralDensityEstimator implements PowerSpectralDensityEstimator {

	/**
	 * Calculates the power spectrum density (PSD) of the given
	 * rr-interval-data. The sampling rate is used to sample the interpolated
	 * rr-interval-data.
	 * 
	 * @param rr
	 *            rr-interval-data in seconds
	 * @return power spectrum density (PSD) of RR-interval-data (unit is s squared)
	 */
	@Override
	public PowerSpectrum calculateEstimate(RRData rr) {
		double[] power = calculatePower(rr);
		double[] frequencies = calculateFrequencies(rr);
		
		return new PowerSpectrum(power, frequencies);
	}
	
	private double[] calculatePower(RRData rr) {
		double[] rrY = rr.getValueAxis();
		var calc = new AvgSampleSizeCalculator();
		double avgSampleSize = calc.process(rr).getValue();
		
		final var fft = new FastFourierTransformer(DftNormalization.STANDARD);
		final Complex[] fftResult = fft.transform(rrY, TransformType.FORWARD);

		// Calculate the power spectral density of the Fourier Transform
		var betrag = new double[rrY.length];
		for (var i = 0; i < rrY.length; i++) {
			//Calculate power spectral density
			//Unit: s*s / Hz = s*s*s
			double abs = fftResult[i].abs();
			betrag[i] =  abs * abs / (rrY.length * (1 / avgSampleSize) * 0.5); 
		}
		
		return betrag;
	}
	
	private double[] calculateFrequencies(RRData rr) {
		var calc = new AvgSampleSizeCalculator();
		double sampleFrequency = 1 / calc.process(rr).getValue();
		double maxSampleFrequency = 0.5 * sampleFrequency;

		// Calculate Frequencies to corresponding power values
		// Frequency Steps in Hz
		double frequencySteps = (2 * maxSampleFrequency) / (rr.getTimeAxis().length);
		var frequencies = new double[rr.getTimeAxis().length];

// Path: lib/src/main/java/hrv/lib/hrv/RRData.java
// public class RRData {
// 
// 	private final TimeUnit rrValueAxisUnit;
// 	private final /*@ spec_public */ TimeUnit rrTimeAxisUnit;
// 
// 	private final double[] rrValuesAxis;
// 	private final double[] rrTimeAxis;
// 
// 	public RRData(double[] timeAxis, TimeUnit timeAxisUnit, double[] valueAxis, TimeUnit valueAxisUnit) {
// 		this.rrValuesAxis = valueAxis;
// 		this.rrValueAxisUnit = valueAxisUnit;
// 
// 		this.rrTimeAxis = timeAxis;
// 		this.rrTimeAxisUnit = timeAxisUnit;
// 	}
// 
// 	/**
// 	 * Creates a new RRData object from an RR-Interval
// 	 * 
// 	 * @param rrInterval
// 	 *            RR-Interval
// 	 * @param unit
// 	 *            of the given RR-Data.
// 	 * @return new RRData object with given RR-Interval data.
// 	 */
// 	public static RRData createFromRRInterval(double[] rrInterval, TimeUnit unit) {
// 		var rrTimeAxis = new double[rrInterval.length];
// 		for (var i = 1; i < rrInterval.length; i++) {
// 			rrTimeAxis[i] = rrTimeAxis[i - 1] + rrInterval[i - 1];
// 		}
// 
// 		return new RRData(rrTimeAxis, unit, rrInterval, unit);
// 	}
// 
// 	public double[] getValueAxis() {
// 		return rrValuesAxis;
// 	}
// 
// 	public double[] getTimeAxis() {
// 		return rrTimeAxis;
// 	}
// 	
// 	public TimeUnit getValueAxisUnit() {
// 		return rrValueAxisUnit;
// 	}
// 
// 	public TimeUnit getTimeAxisUnit() {
// 		return rrTimeAxisUnit;
// 	}
// }
// 
// Path: lib/src/main/java/hrv/lib/hrv/calc/parameter/AvgSampleSizeCalculator.java
// public class AvgSampleSizeCalculator implements HRVDataProcessor {
// 
// 	@Override
// 	public HRVParameter process(RRData data) {
// 		double maxX = data.getTimeAxis()[data.getTimeAxis().length - 1];
// 		double minX = data.getTimeAxis()[0];
// 		double xLength = maxX - minX;
// 		return new HRVParameter(HRVParameterEnum.AVG_SAMPLE_SIZE, xLength / (data.getTimeAxis().length - 1),
// 				"1 / " + data.getTimeAxisUnit().toString());
// 	}
// 	
// 
// 	@Override
// 	public boolean validData(RRData data) {
// 		return data.getValueAxis().length > 0;
// 	}
// }


// Path: lib/src/main/java/hrv/lib/hrv/calc/psd/StandardPowerSpectralDensityEstimator.java
import hrv.lib.hrv.RRData;
import hrv.lib.hrv.calc.parameter.AvgSampleSizeCalculator;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.transform.DftNormalization;
import org.apache.commons.math3.transform.FastFourierTransformer;
import org.apache.commons.math3.transform.TransformType;

package hrv.lib.hrv.calc.psd;


/**
 * Uses FFT to estimate the power spectra density of RR-Data-Intervals
 * 
 * @author Julian
 *
 */
public class StandardPowerSpectralDensityEstimator implements PowerSpectralDensityEstimator {

	/**
	 * Calculates the power spectrum density (PSD) of the given
	 * rr-interval-data. The sampling rate is used to sample the interpolated
	 * rr-interval-data.
	 * 
	 * @param rr
	 *            rr-interval-data in seconds
	 * @return power spectrum density (PSD) of RR-interval-data (unit is s squared)
	 */
	@Override
	public PowerSpectrum calculateEstimate(RRData rr) {
		double[] power = calculatePower(rr);
		double[] frequencies = calculateFrequencies(rr);
		
		return new PowerSpectrum(power, frequencies);
	}
	
	private double[] calculatePower(RRData rr) {
		double[] rrY = rr.getValueAxis();
		var calc = new AvgSampleSizeCalculator();
		double avgSampleSize = calc.process(rr).getValue();
		
		final var fft = new FastFourierTransformer(DftNormalization.STANDARD);
		final Complex[] fftResult = fft.transform(rrY, TransformType.FORWARD);

		// Calculate the power spectral density of the Fourier Transform
		var betrag = new double[rrY.length];
		for (var i = 0; i < rrY.length; i++) {
			//Calculate power spectral density
			//Unit: s*s / Hz = s*s*s
			double abs = fftResult[i].abs();
			betrag[i] =  abs * abs / (rrY.length * (1 / avgSampleSize) * 0.5); 
		}
		
		return betrag;
	}
	
	private double[] calculateFrequencies(RRData rr) {
		var calc = new AvgSampleSizeCalculator();
		double sampleFrequency = 1 / calc.process(rr).getValue();
		double maxSampleFrequency = 0.5 * sampleFrequency;

		// Calculate Frequencies to corresponding power values
		// Frequency Steps in Hz
		double frequencySteps = (2 * maxSampleFrequency) / (rr.getTimeAxis().length);
		var frequencies = new double[rr.getTimeAxis().length];

		for (var i = 0; i < frequencies.length; i++) {