/* -*- c++ -*- */
/* 
 * Copyright 2015 <+YOU OR YOUR COMPANY+>.
 * 
 * This is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3, or (at your option)
 * any later version.
 * 
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this software; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <gnuradio/io_signature.h>
#include "howto_detect_ff_impl.h"
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>

namespace gr {
  namespace howto {

    howto_detect_ff_sptr
    howto_make_detect_ff(float pfa, int L, int samples)
    {
      return gnuradio::get_initial_sptr
        (new howto_detect_ff_impl(pfa, L, samples));
    }

    static const int MIN_IN = 1; // mininum number of input streams
    static const int MAX_IN = 1; // maximum number of input streams
    static const int MIN_OUT = 1; // minimum number of output streams
    static const int MAX_OUT = 1; // maximum number of output streams
    float mem = 0; //Global Variable
     
     /*
     * The private constructor
     */
    
    
    howto_detect_ff_impl::howto_detect_ff_impl (float pfa, int L, int samples)
      : gr::block("howto_detect_ff",
              gr::io_signature::make (MIN_IN, MAX_IN, sizeof (float)),
              gr::io_signature::make (MIN_OUT, MAX_OUT, sizeof(float)),
       d_pfa(pfa), d_L(L), d_samples(samples)
    {
    }

    /*
     * Our virtual destructor.
     */
    howto_detect_ff_impl::~howto_detect_ff_impl()
    {
    }

//------------functions---------------

/*This function gives the CDF value for the pfa in input*/

float TracyWidom (float p){

float pd, tw;
tw = 0.45;
pd = 1 - p;
if (pd >= 0.01 && pd <= 0.05){
tw = 18*(pd - (17/75)); printf("a - %f\n", tw);
}else if (pd >= 0.05 && pd <= 0.1){
tw = 8*(pd - (179/400)); printf("b - %f\n", tw);
}else if (pd >= 0.1 && pd <= 0.3){
tw = (87/20)*(pd - (643/870)); printf("c - %f\n", tw);
}else if (pd >= 0.3 && pd <= 0.5){
tw = (16/5)*(pd - (287/320)); printf("d - %f\n", tw);
}else if (pd >= 0.5 && pd <= 0.7){
tw = (17/5)*(pd - (297/340)); printf("e - %f\n", tw);
}else if (pd >= 0.7 && pd <= 0.9){
tw = (5.2)*(pd - (0.813)); printf("f - %f\n", tw);
}else if (pd >= 0.9 && pd <= 0.95){
tw = (53/5)*(pd - (909/1060)); printf("g - %f\n", tw);
}else if (pd >= 0.95 && pd <= 1){
tw = 26*(pd - (593/650)); printf("h - %f\n", tw);
}else{
printf ("wrong pfa value: it must be between 0 and 1\n");
printf ("the pfa value standard is 0.1\n");
tw = (53/5)*(0.9 - (909/1060));
}
return tw;
}
/*this function calculates the threshold for the test
the inputs are: ns = numbers of samples;
L = length of the correlation function;
pfa = probability of false alarm*/
float gamma (int ns, int L, float pfa){
  float A, B, C, ratio1, ratio2, g;
  A = sqrt(ns)+sqrt(L);
  B = sqrt(ns)-sqrt(L);
  C = ns*L;
  ratio1 = pow(A,2)/pow(B,2);
  ratio2 = 1+( pow(A,-0.667) / pow(C,0.167) )*pfa;
  g = ratio1*ratio2;

  return g;
}
/*This function makes the detection test*/

float test (float r, float t){
    float decision;
    if(r > -1){
      if(r <= t){
        decision = 0;
       }else{
        decision = 1;
}
}
return decision;}
//-------------end functions-----------

    
    int
    howto_detect_ff_impl::general_work (int noutput_items,
                       gr_vector_int &ninput_items,
                       gr_vector_const_void_star &input_items,
                       gr_vector_void_star &output_items)
    {
        const float *in = (const float *) input_items[0];
        float *out = (float *) output_items[0];
        float vett [d_samples];
        int lenght = floor(d_samples / d_L) * d_L;
        int p, j, k;
        float thr, lmax, lmin, ratio, TW;
        // char c[1];
        gsl_matrix * hankel = gsl_matrix_alloc (lenght,d_L);
        gsl_matrix * V = gsl_matrix_alloc (d_L,d_L);
        gsl_vector * S = gsl_vector_alloc (d_L);
        gsl_vector * temp = gsl_vector_alloc (d_ L);
        FILE *story;
        k=0; ratio = -1;
        gsl_matrix_set_zero (hankel);
        TW = TracyWidom (d_pfa);
        thr = gamma(lenght, d_L, TW);
	for (int i = 0; i < noutput_items; i++){
         vett[k]=in[i];
         k++;

         if (k >= lenght){
            k = 0;
            story = fopen("filestory.txt", "a");
            for( p=0; p<lenght; p++ ){
                for( j=0; j<d_L; j++ ){
                    gsl_matrix_set (hankel, p, j, vett[p+j]);
}
}
                gsl_linalg_SV_decomp (hankel, V, S, temp);
                lmax = gsl_vector_get(S, 0);
                lmin = gsl_vector_get(S, (d_L -1));
                ratio = lmax/lmin;
                mem = test(ratio, thr);
                fprintf(story, "%f - ratio=%f - soglia=%f\n ", mem, ratio, thr);
                fclose(story);
}
                out[i] = mem;
}
        // Do <+signal processing+>
        // Tell runtime system how many input items we consumed on
        // each input stream.
        consume_each (noutput_items);

        // Tell runtime system how many output items we produced.
        return noutput_items;
    }

}   /* namespace howto */
} /* namespace gr */