package com.diampark.test.Class;

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;

public class calcul_PSD {

	public final static FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);

	public static double[] calc_PSD_non_normalise(double[] tab_valeurs) {

		Complex[] fft_valeurs = fft.transform(tab_valeurs,TransformType.FORWARD);

		int n = fft_valeurs.length;
		n/=2;

		double[] PSD_valeurs = new double[n];
		for (int i=0; i<n; i++) PSD_valeurs[i]=fft_valeurs[i].abs()*fft_valeurs[i].abs();

		return PSD_valeurs;
	}

	public static boolean is_power_of_two(int n) {
		if (n==1) return true;

		else return (n%2==0)&(is_power_of_two(n/2));
	}

	public static double[] ajuster_taille(double[] tab_valeurs) {
		int n = tab_valeurs.length;
		int p =n;

		while (!is_power_of_two(p)) {
			p++;
		}

		double[] tab_2 = new double[p];

		for (int i=0;i<n;i++) tab_2[i]=tab_valeurs[i];
		for (int i=n;i<p;i++)tab_2[i]=0.0;

		return tab_2;
	}

	public static double[] calcul_puissance_1(double[] tab_valeurs,double freq_echantillonnage,double[] fenetre) {

		double[] PSD_valeurs = calc_PSD_non_normalise(ajuster_taille(tab_valeurs));
		int n = PSD_valeurs.length;

		double[] freqs = new double[n];
		for (int i=0;i<n;i++) freqs[i]=i*freq_echantillonnage/(2*n);

		double pas = freq_echantillonnage/(2*n);

		double integr_tot = 0.0;
		double integr_fenetre=0.0;

		for (int i=0;i<n-1;i++) {
			integr_tot+=pas*(PSD_valeurs[i+1]+PSD_valeurs[i])/2.0;

			if (freqs[i]>=fenetre[0]&freqs[i+1]<=fenetre[1]) integr_fenetre+=pas*(PSD_valeurs[i+1]+PSD_valeurs[i])/2.0;
		}

		return new double[]{integr_fenetre,integr_tot};
	}

	public static double[] calcul_puissance_2(double[] tab_valeurs,double freq_echantillonnage,double[] fenetre) {

		double[] PSD_valeurs = calc_PSD_non_normalise(ajuster_taille(tab_valeurs));
		int n = PSD_valeurs.length;

		double[] freqs = new double[n];
		for (int i=0;i<n;i++) freqs[i]=i*freq_echantillonnage/(2*n);

		for (int i=0;i<n;i++) {
			System.out.print(freqs[i]+",");
		}
		System.out.println("");
		for (int i=0;i<n;i++) {
			System.out.print(PSD_valeurs[i]+",");
		}
		System.out.println("");

		double pas = freq_echantillonnage/(2*n);

		double pmax=0.0;
		double freqmax=0.0;

		for (int i=0;i<n;i++) {
			if (freqs[i]>=fenetre[0]&freqs[i]<=fenetre[1]) {
				if (PSD_valeurs[i]>pmax) {
					pmax=PSD_valeurs[i];
					freqmax=freqs[i];
				}
			}
		}

		double integr_pic = 0.0;
		double integr_fenetre=0.0;

		for (int i=0;i<n-1;i++) {

			if (freqs[i]>=fenetre[0]&freqs[i+1]<=fenetre[1]) {
				integr_fenetre+=pas*(PSD_valeurs[i+1]+PSD_valeurs[i])/2.0;
				if ((freqs[i]>=freqmax-1)&(freqs[i+1]<=freqmax+1)) integr_pic+=pas*(PSD_valeurs[i+1]+PSD_valeurs[i])/2.0;
			}
		}

		return new double[] {integr_pic,integr_fenetre};
	}

	public static double calcul_Q1(double[] tab1,double[] tab2,double[] tab3,double freq_echantillonnage,double[] fenetre) {

		double[] p1 = calcul_puissance_1(tab1,freq_echantillonnage,fenetre);
		double[] p2 = calcul_puissance_1(tab2,freq_echantillonnage,fenetre);
		double[] p3 = calcul_puissance_1(tab3,freq_echantillonnage,fenetre);

		return (p1[0]+p2[0]+p3[0])/(p1[1]+p2[1]+p3[1]);
	}

	public static double calcul_Q2(double[] tab1,double[] tab2,double[] tab3,double freq_echantillonnage,double[] fenetre) {

		double[] p1 = calcul_puissance_2(tab1,freq_echantillonnage,fenetre);
		double[] p2 = calcul_puissance_2(tab2,freq_echantillonnage,fenetre);
		double[] p3 = calcul_puissance_2(tab3,freq_echantillonnage,fenetre);

		return (p1[0]+p2[0]+p3[0])/(p1[1]+p2[1]+p3[1]);
	}

	public static double[] calcul_Qs(double[] tab1,double[] tab2,double[] tab3,double freq_echantillonnage,double[] fenetre) {
		return new double[] {calcul_Q1(tab1,tab2,tab3,freq_echantillonnage,fenetre),calcul_Q2(tab1,tab2,tab3,freq_echantillonnage,fenetre)};
	}

}