#include "dvbt_signauxtests_v02.h"

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>

/*!	\fn	CreateDVBTSignauxtests()

	\brief
		\~english	Memory allocation and initialisation of structure DVBTSignauxtests.
		\~french	Allocation de mémoire et initialisation de la structure DVBTSignauxtests.

	\return
		\~english	Pointer on structure DVBTSignauxtests.
		\~french	Pointeur sur structure DVBTSignauxtests.
	   \~french   A liberer. NULL = erreur d'allocation.

*/
DVBTSignauxtests *CreateDVBTSignauxtests (void)
{
   /* Memory allocation for structure DVBTSignauxtests. */
   DVBTSignauxtests *DS = malloc (sizeof *DS);
   if (DS == NULL)
   {
      printf
         ("\nERROR !\tCreateDVBTSignauxtests : unable to allocate memory for DVBTSignauxtests structure.\n\n");
      exit (EXIT_FAILURE);
   }
   else
   {
      /* all elements to zero */
      static const DVBTSignauxtests z = { 0, 0, 0, 0, 0 };
      *DS = z;
   }
   return DS;
}

/*!	\fn	FreeDVBTSignauxtests(DVBTSignauxtests *DS)

	\brief
		\~english	Memory release of structure DVBTSignauxtests.
		\~french	Libération de la mémoire de la structure DVBTSignauxtests.

	\param	DS :
		\~english	Pointer on structure DVBTSignauxtests.
		\~french	Pointeur sur structure DVBTSignauxtests.

	\return
		-.
*/
void FreeDVBTSignauxtests (DVBTSignauxtests * DS)
{
   /* Test structure existence. */
   if (DS == NULL)
   {
      printf
         ("\nWARNING !\tFreeDVBTSignauxtests : structure DVBTSignauxtests does not exist.\n\n");
      return;
   }

   /* Memory release. */
   free (DS);
}

/*!	\fn	UseDVBTSignauxtests(DVBTSignauxtests *DS, int source_reelle[],int source_imaginaire[],int surechantillonage,int echantillon)

	\brief
		\~english	Memory allocation and initialisation of structure DVBTSignauxtests
		\~french	Allocation mémoire et initialisation de la structure DVBTSignauxtests

	\param	DS                   :  Pointeur sur structure DVBTSignauxtests.
	\param	source_reelle[]      :  Tableau correspondant à une séquence (+1 -1) correspondant à une partie reelle
	\param	source_imaginaire[]  :  Tableau correspondant à une séquence correspondant à une partie imaginaire
	\param	surechantillon       :  Entier correspondant à un facteur de suréchantillonage
	\param	echantillon          :  Entier correspondant au nombre d'échantillons du signal source
	\return
		-.

*/
void UseDVBTSignauxtests (DVBTSignauxtests * DS, int const source_reelle[],
                          int const source_imaginaire[],
                          int surechantillonage, int echantillon)
{
   /* récupération du facteur de suréchantillonnage */
   int const Ns = surechantillonage;
   /* récupération du nombre d'échantillons du signal source */
   int const N = echantillon;
   double *upsignal_re = NULL;
   double *upsignal_im = NULL;

/*------------------------------------------------------------------------------------------------*/
/*                                Surechantillonnage des signaux                                  */
/*------------------------------------------------------------------------------------------------*/

/*Suréchantillonnage du signal partie réelle par un facteur Ns qui a été demandé à l'utilisateur*/
   upsignal_re = UpsamplingDVBTSignauxtests (source_reelle, Ns, N);
   printf ("\nSignal partie reelle surechantillonnée de %d :", Ns);
/*Affichage du signal partie réelle surechantillonné*/
   {
      int i;

      for (i = 0; i < N * Ns; i++) /*La taille du signal surechantillonnée sera le résultat de la multiplication */
      {                         /*du nombre d'échantillons au départ par le facteur de surechantillonnage: N*Ns */
         printf ("%f ", upsignal_re[i]);
      }
      printf ("\n");
   }
/*Suréchantillonnage du signal partie imaginaire*/
   upsignal_im = UpsamplingDVBTSignauxtests (source_imaginaire, Ns, N);
   printf ("\nSignal partie imaginaire surechantillonnée de %d :", Ns);
/*Affichage du signal partie imaginaire surechantillonné*/

   {
      int i;
      for (i = 0; i < N * Ns; i++) /*La taille du signal surechantillonnée sera le résultat de la multiplication */
      {                         /*du nombre d'échantillons au départ par le facteur de surechantillonnage: N*Ns */
         printf ("%f ", upsignal_im[i]);
      }
      printf ("\n");
   }
/*------------------------------------------------------------------------------------------------*/
/*                Calcul des coefficients du filtre à racine carrée de Nyquist                    */
/*------------------------------------------------------------------------------------------------*/
   {
      double periodes;
      double points;
      double roll_off;
      printf
         ("\n\nCalcul des coefficients du filtre à racine carrée de Nyquist :");
      printf ("\nNombre de périodes :");
      scanf ("%lf", &periodes);

      printf ("\nNombre de points par périodes :");
      scanf ("%lf", &points);

      printf ("\nroll Off:");
      scanf ("%lf", &roll_off);
      {
         double *signal_re_conv_nyquist = NULL;
         double *signal_im_conv_nyquist = NULL;
         {
            double *coef_nyquist =
               RaisedCosineSQRT (periodes, points, roll_off);

/*-------------------------------------------------------------------------------------------------*/
/*Convolution des signaux surechantillonnés et les coefficients du filtre à racine carré de Nyquist*/
/*-------------------------------------------------------------------------------------------------*/
/*filtre : coef_nyquist , taille = periodes*points*/
/*signal : upsignal (signal surechantillonné) , taille = N*Ns (echantillon_source*facteur_surechantillonage)*/

/*Signal partie reelle*/
            signal_re_conv_nyquist =
               conv (coef_nyquist, upsignal_re, N * (Ns), periodes * points);
            if (signal_re_conv_nyquist != NULL)
            {
               printf
                  ("\nSignal partie reelle sortie du filtre de Nyquist:\n");
/*Affichage du signal partie reelle après le passage dans le filtre à racine de Nyquist*/
               {
                  int i;
                  for (i = 0; i < (periodes * points + N * Ns - 1); i++)
                  {
                     printf ("%f ", signal_re_conv_nyquist[i]);
                  }
               }
/*Signal partie imaginaire*/
               signal_im_conv_nyquist =
                  conv (coef_nyquist, upsignal_im, N * (Ns),
                        periodes * points);
            }
            else
            {
               printf ("signal_re_conv_nyquist = NULL\n");
            }
         }

/*------------------------------------------------------------------------------------------------*/
/*                                       Filtre FIR                                               */
/*------------------------------------------------------------------------------------------------*/
         printf ("\n\nFiltre FIR");

/*------------------------------------------------------------------------------------------------*/
/*                                    CIC interpolateur                                           */
/*------------------------------------------------------------------------------------------------*/
         {
            int nombre_etages;
            int interp;

            printf ("\nEntrez le Nombre d'étages du CIC :");
            scanf ("%d", &nombre_etages);
            printf ("\nEntrez le Rapport d'interpolation :");
            scanf ("%d", &interp);
            {
               int taille_signal = (periodes * points + N * Ns - 1) * interp;

               //  double sortie_interp[taille_cic_interpolateur];

               double *sortie_interp =
                  malloc (sizeof *sortie_interp * taille_signal);

               CicNum_stages (DS, nombre_etages); /*Initialisation de la structure du CIC en fonction du nombre d'étages */
               printf ("\n");

               printf ("\navant interpolation :\n");
               {
                  int i;
                  for (i = 0; i < (periodes * points + N * Ns - 1); i++)
                  {
                     printf (" %f", signal_re_conv_nyquist[i]);
                  }
               }
               printf ("\n");
               {
                  int i;
                  for (i = 0; i < taille_signal; ++i)
                  {
                     double interpolation = CicInterpolate (DS, signal_re_conv_nyquist[i], 1); /*Utilisation du CIC Interpolation */
                     sortie_interp[i * interp] = interpolation;
                     int j;
                     for (j = 1; j < interp; j++)
                     {
                        sortie_interp[j + interp * i] = interpolation;
                     }
                  }
               }
               printf ("\ninterpolation :\n");
               {
                  int i;
                  for (i = 0; i < taille_signal; i++)
                  {
                     printf (" %lf", sortie_interp[i]);
                  }
               }

/*-------------------------------------------------------------------------------------------------*/
/*                                Transposition fréquentielle                                      */
/*-------------------------------------------------------------------------------------------------*/
               printf ("\nappel fonction transposition");
               {
                  double pas_echantillon = 1 / ((double) (Ns * interp));

                  FrequencyTransposition (sortie_interp, pas_echantillon,
                                          taille_signal);
               }
            }
         }
      }
   }
}

/*!	\fn	double * UpsamplingDVBTSignauxtests(double sig[],int fact_surechantillon,int echantillon)
	\brief
		\~english	.
		\~french	  Fonction permettant de suréchantillonner un signal

	\param	sig                  :  Signal que l'on veut échantillonner.
	\param	fact_surechantillon  :  Facteur de suréchantillonnage
	\param	echantillon          :  Nombre d'échantillons avant le suréchantillonnage

	\return
	   \~french   Pointeur sur tableau du signal surechantillonné.
	   \~french   A liberer. NULL = erreur d'allocation.

*/

double *UpsamplingDVBTSignauxtests (int const sig[],
                                    int fact_surechantillon, int echantillon)
{
   /* récupération du facteur de suréchantillonnage */
   int const Ns = fact_surechantillon;
   int const N = echantillon;
   double *up_signal = malloc (sizeof *up_signal * (N * Ns));
   if (up_signal != NULL)
   {
/*-----------------surechantillonage----------------------*/
      /*pointeur pour parcourir le tableau correspndant au signal */
      int const *p;
      int i;
      for (i = 0, p = sig; i < N * Ns; ++i)
      {
         /*A chaque indice non multiple du facteur d'échantillonnage+1 */
         if (i % Ns)
         {
            /*on insère les valeurs 0 */
            up_signal[i] = 0;
         }
         else
         {
            /*sinon on insère les valeurs du signal originales */
            up_signal[i] = *p;
            p++;
         }
      }
   }
   return up_signal;
}

/*!	\fn	double * conv(double *filtre, double *signal, int taille_entree, int taille_filtre)

	\brief
		\~english	Convolution of x[N] with H[N]
		\~french	Fonction réalisant la convolution de X[N] avec H[N], le résultat sera stocké dans Y[N]

	\param	filtre               :	tableau correspondant au filtre.
	\param	signal               :	tableau correspondant au signal d'entree.
	\param	taille_entree        :  taille du signal d'entree.
    \param	taille_filtre        :  taille du filtre.

	\return
      \~french   Pointeur sur tableau correspondant au signal convolué.
	   \~french   A liberer. NULL = erreur d'allocation.

*/
double *conv (double const filtre[], double const signal[], int taille_entree,
              int taille_filtre)
{
   double *y = NULL;
   /*Récupération de la taille du signal d'entree */
   int const L = taille_entree;
   /*Récupération de la taille du filtre */
   int const M = taille_filtre;
   /*Déclaration de tableau de taille correspondant à la longueur de la
      sortie du signal convoluée */
   double *x = malloc (sizeof *x * (L + M));

   printf ("conv(%p, %p, %d, %d)\n", filtre, signal, taille_entree,
           taille_filtre);

   if (x != NULL)
   {
      double *h = malloc (sizeof *h * (L + M));
      {
         int const nby = L + M - 1;
         /*Allocation en mémoire */
         y = malloc (sizeof *y * nby);

         if (y != NULL)
         {
            int i;
            for (i = 0; i < nby; i++)
            {
               x[i] = 0;        /*Rajout de 0 dans le vecteur signal d'entree */
               h[i] = 0;        /*Rajout de 0 dans le vecteur filtre */
            }

            for (i = 0; i < M; i++)
            {
               h[i] = filtre[i]; /*Récupération des coefficcients du filtre */
            }

            for (i = 0; i < L; i++)
            {
               x[i] = signal[i]; /*Récupération du signal d'entree */
            }

            for (i = 0; i < nby; i++) /*Initialisation du tableau de sortie */
            {
               y[i] = 0;
            }

            {
               int n;

               for (n = 0; n < nby; n++) /*Calcul de la convolution entre les 2 signaux */
               {
                  int m;
                  for (m = 0; m <= n; m++)
                  {
                     assert (n < L + M);
                     assert (m < L + M);
                     assert (n - m < L + M);
                     assert (n - m >= 0);

                     y[n] = h[m] * x[n - m] + y[n];
                  }
               }
            }
         }
         free (h), h = 0;
      }
      free (x), x = 0;
   }
   else
   {
      printf ("conv() : x=NULL\n");

   }
   return y;
}

/*!	\fn	void FrequencyTransposition(double *signal_reel, double *signal_imag, double pas_echantillon, int taille_signal)

	\brief
		\~english
		\~french	Fonction réalisant une transposition fréquentielle

	\param	signal_reel               :	tableau correspondant au signal partie reelle.
	\param	signal_imag               :	tableau correspondant au signal partie imaginaire.
	\param	pas_echantillon           : pas d'échantillons.
	\param	taille_signal             : taille des tableaux.

   	\return
		-.

*/
void FrequencyTransposition (double const signal_reel[], double tk,
                             int taille_signal)
{
   double fp;

   printf ("\n\nFréquence de la porteuse en MHz :");
   scanf ("%lf", &fp);

   printf
      ("\n\nSignal partie réelle transposé par une porteuse de %f MHz :", fp);

   {
      FILE *FicSortie = fopen ("reel.txt", "w");
      if (FicSortie != NULL)
      {
         int const N = taille_signal; /*récupération de la taille des tableaux */
         int i;
         for (i = 0; i < N; ++i)
         {
            /*Multiplication du signal reel par cos(2*pi*fp*tk) */
            double cosinus = cos (2 * M_PI * fp * i * tk);
            double signal_re_transpose = signal_reel[i] * cosinus;
            fprintf (FicSortie, "%f\n ", signal_re_transpose);
         }
         fclose (FicSortie);
      }
   }
   return;
}

/*!	\fn	double * RaisedCosineSQRT(N,D,e)
	\brief
		\~english
		\~french	  Création d'un vecteur de reel de forme racine d'un cosinus sureleve

	\param	N  : Nombre de périodes.
	\param	D  : Nombre de points par périodes.
	\param	e  : Valeur de Roll-off.

	\return
		\~french   Pointeur sur tableau correspondant aux coefficients réels de forme surcos.
	   \~french   A liberer. NULL = erreur d'allocation.

*/

double *RaisedCosineSQRT (double nb_periodes, double nb_points,
                          double roll_off)
{

/*------déclaration et initialisation des variables------*/
   double const N = nb_periodes;
   double const D = nb_points;
   int const nbw = nb_periodes * nb_points;

   double *w = malloc (sizeof *w * nbw); /*Allocation en mémoire */
   if (w != NULL)
   {
      double const e = roll_off;
      int n;

      double racineD = sqrt (D);

      for (n = 0; n < nbw; n++)
      {
         w[n] = 0;              /*vecteur de 0 */
      }

      for (n = 0; n < nbw; n++)
      {
         int const t = n - N * D / 2;
         double const x = t / D;
         double const oneplus = (1.0 + e) * M_PI / D;
         double const oneminus = (1.0 - e) * M_PI / D;

         if (t == 0.000000)
         {
            w[n] = (4 * e / M_PI + 1 - e) / racineD;
         }
         else if (e == 0.000000)
         {
            w[n] = racineD * sin (M_PI * x) / (3.1416 * t);
         }
         else
         {
            double const den = t * t * 16 * e * e - D * D;

            if (den == 0)
            {
               w[n] =
                  (D * racineD / (8 * e * 3.1416 * t)) * (oneplus *
                                                          sin (oneplus * t) -
                                                          (oneminus * D /
                                                           (4 * e * t)) *
                                                          cos (oneminus * t) +
                                                          (D /
                                                           (4 * e * t * t)) *
                                                          sin (oneminus * t));
            }
            w[n] =
               (4 * e / (3.1416 * racineD)) * (cos (oneplus * t) +
                                               sin (oneminus * t) / (x * 4 *
                                                                     e)) /
               (1.0 - 16 * e * e * x * x);
         }
      }

      for (n = 0; n < N * D; n++)
      {
         printf ("\nw[i]%f", w[n]);
      }
   }
   return w;
}

/*!	\fn	void  CicNum_stages(DVBTSigauxtests *DS, int n)
	\brief
		\~english
		\~french	  Fonction permettant d'initialiser le nombre d'étages du filtre CIC

	\param	DC      : Pointeur sur structure DVBTCic
	\param	n       : Entier correspondant au nombre d'étages du filtre CIC

	\return
		\~french   .

*/

void CicNum_stages (DVBTSignauxtests * DS, int n)
{

   DS->stages = n;

   free (DS->nacc);
   free (DS->diff);
   free (DS->prev);

   /*Allocation mémoire correspondant au nombre d'étages du CIC */
   DS->nacc = malloc (sizeof *DS->nacc * n);
   DS->diff = malloc (sizeof *DS->diff * n);
   DS->prev = malloc (sizeof *DS->prev * n);

   return;
}

/*!	\fn	double  CicInterpolate(DVBTCic *DC, double in, int dump)
	\brief
		\~english
		\~french	  Fonction permettant d'effectuer une interpolation

	\param	DC      : Pointeur sur structure DVBTCic
	\param	in      : Double représentant les échantillons en entrées avant interpolation.
	\param	dump    : Entier représentant le taux d'interpoaltion.

	\return
		\~french     Double représentant les échantillons après interpolation

*/
double CicInterpolate (DVBTSignauxtests * DS, double in, int dump)
{
   int i;

   /*Dump autorise l'entrée d'un nouvel échantillon */
   if (dump)
   {
      DS->diff[0] = in - DS->prev[0];
      DS->prev[0] = in;
      for (i = 1; i < DS->stages; i++)
      {
         DS->diff[i] = DS->diff[i - 1] - DS->prev[i];
         DS->prev[i] = DS->diff[i - 1];
      }
      DS->nacc[0] += DS->diff[DS->stages - 1];
   }

   for (i = 0; i < (DS->stages - 1); i++)
   {
      DS->nacc[i + 1] += DS->nacc[i];
   }

   return (DS->nacc[DS->stages - 1]);

}

/*!	\fn	ResetDVBTSignauxtests(DVBTSignauxtests *DS)

	\brief
		\~english	Reset of structure DVBTSignauxtests.
		\~french	Ré-initialisation de la structure DVBTSignauxtests.

	\param	DS :
		\~english	Pointer on structure DVBTSignauxtests.
		\~french	Pointeur sur structure DVBTSignauxtests.

	\return
		-.

*/
void ResetDVBTSignauxtests (DVBTSignauxtests * DS)
{

   return;
}