/* See dm_spectrum.h */
/* Last edited on 2008-04-19 11:28:48 by stolfi */

#define dm_spectrum_C_COPYRIGHT \
  "Copyright © 2005  by the State University of Campinas (UNICAMP)" \
  " and the Federal Fluminense University (UFF)"

#include <dm_spectrum.h>
#include <dm_nucleic.h>
#include <dm_basic.h>
#include <dm_sample.h>
#include <dm_datum.h>
#include <dm_seq.h>

#include <msm_ps_tools.h>
#include <msm_image.h>

#include <vec.h>
#include <nget.h>
#include <fget.h>
#include <filefmt.h>
#include <jsrandom.h>
#include <pswr.h>

#include <fftw3.h>

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

/* IMPLEMENTATIONS */

double_vec_t dm_spectrum_from_seq(dm_seq_t *seqp)
  {
    /* Get sequence size: */
    int nc = dm_CHANNELS;
    int nd = dm_seq_num_datums(seqp); /* Number of sample datums. */

    /* Get scaling factors: */
    double *sfac = seqp->sfac.f; /* Sample scaling factor per channel. */
    
    /* Allocate vectors {in,ot} for {fftw3.h} routines: */
    fftw_complex *in = fftw_malloc(sizeof(fftw_complex)*nd); /* Input vector. */
    fftw_complex *ot = fftw_malloc(sizeof(fftw_complex)*nd); /* Output vector. */

    /* Precompute parameters for {fftw3.h} routines: */
    fftw_plan plan = fftw_plan_dft_1d(nd, in, ot, FFTW_FORWARD, FFTW_ESTIMATE);

    /* Compute the total power spectrum {pwr[0..fMax]}: */
    int fMax = nd/2; /* Maximum frequency in power spectrum. */
    int nf = fMax+1; /* Number of enries in power spectrum. */

    /* Allocate and clear the power spectrum: */
    double_vec_t pwr = double_vec_new(nf);
    int f;
    for (f = 0; f <= fMax; f++) { pwr.e[f] = 0; }
    
    int c;
    for (c = 0; c < nc; c++)
      { /* Copy channel {c} samples, scaled, to {in} vector, as complex numbers: */
        int i;
        for (i = 0; i < nd; i++)
          { double s = dm_seq_get_sample(seqp, i, c);
            in[i][0] = dm_sample_decode(s)*sfac[c];
            in[i][1] = 0; 
          }

        /* Compute Fourier transform: */
        fftw_execute(plan);

        /* Compute the power spectrum and add it to {pwr[0..nf-1]}: */
        double norm = 1/((double)nd); /* Power normalization factor. */
        int f;
        for (f = 0; f <= fMax; f++)
          { /* Add square modulus of Fourier terms {f} and {-f}, if distinct: */
            double P = ot[f][0]*ot[f][0] + ot[f][1]*ot[f][1];
            int g = (nd - f) % nd;
            if (g != f) { P += ot[g][0]*ot[g][0] + ot[g][1]*ot[g][1]; }
            pwr.e[f] += norm*P;
          }
      }
    return pwr;
  }

void dm_spectrum_postscript_plot(msm_ps_tools_t *dp, dm_seq_t *seqp)
  { 
    /* Get sequence size and circularity: */
    int nd = dm_seq_num_datums(seqp); /* Number of sample datums. */
    if (nd == 0) { return; }

    /* Get scaling factors: */
    double *sfac = seqp->sfac.f; /* Sample scaling factor per channel. */
    
    /* Estimate the maximum power {pwrMax} per specral frequency: */
    int c;
    double smpUMax = dm_sample_decode(dm_sample_VALID_MAX);  /* Max value of an unscaled sample. */
    double datPwrMax = 0.0;      /* Max power per datum. */
    for (c = 0; c < dm_CHANNELS; c++)
      { double smpMax = smpUMax * sfac[c]; datPwrMax += smpMax*smpMax; }
    /* Assume that power spectrum is more or less uniform: */
    double pwrMax = 2*datPwrMax;

    /* Choose the nominal Y range {[yMin_yMax]}: */
    double yMin = 0.0;
    double yMax = 1.50 * pwrMax;
    
    /* Compute the sequence's power spectrum: */
    double_vec_t pwr = dm_spectrum_from_seq(seqp);
      
    /* Plot the graphs: */
    msm_ps_tools_draw_histogram(dp, pwr.ne, NULL, pwr.e, yMin, yMax); 
    
    free(pwr.e);
  }

void dm_spectrum_postscript_plot_named
  ( dm_seq_t *seqp, 
    double hSize, 
    double vSize,
    double fontSize,
    char *name,
    char *tag
  )
  { /* Create the {msm_ps_tools_t} object stream from it: */
    msm_ps_tools_t *dp = msm_ps_tools_new_graph
      ( NULL, name, tag,
        hSize, vSize, 
        /*scaleL*/  TRUE,  /*titleL*/ FALSE,
        /*scaleR*/  FALSE, /*titleR*/ FALSE,
        /*scaleB*/  TRUE,  /*titleB*/ FALSE,
        /*scaleT*/  FALSE, /*titleT*/ FALSE,
        fontSize,
        /*maxLabChars*/ 7,
        msm_EPS_MARGIN_MM
      );
    
    /* Plot the sequence: */
    dm_spectrum_postscript_plot(dp, seqp);
    /* Close and cleanup: */
    msm_ps_tools_close(dp);
  }