/*  
 * spectrum_analyzer.c -- Analizador de espectros.  
 *  
 * Este fichero fuente puede encontrarse en:  
 * http://www.ace.ual.es/~vruiz/docencia/redes/practicas/spectrum_analyzer.c  
 *  
 * Compilar escribiendo (el paquete fftw debería estar instalado!):  
 * gcc spectrum_analyzer.c -o spectrum_analyzer spin.o -lfftw3 -lm  
 *  
 * Más información en: http/www.fftw.org  
 *  
 * gse. 2007  
 */  
 
#include <stdio.h>  
#include <stdlib.h>  
#include <math.h>  
#include <fftw3.h>  
#include "spin.h"  
 
main(int argc, char *argv[]) {  
  if(argc < 2) {  
    fprintf(stderr,"%s input_signal.float > spectrum.txt\n",argv[0]);  
  } else {  
    FILE *input_file = fopen(argv[1],"rb");  
    if(!input_file) {  
      fprintf(stderr,"%s: unable to open input file \"%s\"\n"  
              ,argv[0],argv[2]);  
      exit(1);  
    }  
 
    int samples = compute_number_of_samples(input_file);  
    fprintf(stderr,"%s: number of samples = %d\n",argv[0],samples);  
 
    double       *in  = (double       *)fftw_malloc(samples*sizeof(double      ));  
    fftw_complex *out = (fftw_complex *)fftw_malloc(samples*sizeof(fftw_complex));  
 
    /* Creamos un plan para la fftw */  
    fftw_plan p = fftw_plan_dft_r2c_1d(samples, in, out, 0);  
 
    read_signal(in, input_file, samples);  
 
    /* Calculamos la transformada de Fourier */  
    fftw_execute(p);  
 
    /* Generamos el espectro de la señal */ {  
      int k;  
      for(k=0; k<samples/2; k++) {  
double power_spectrum = sqrt(out[k][0]*out[k][0] + out[k][1]*out[k][1]);  
fprintf(stdout,"%f\n",power_spectrum);  
      }  
    }  
 
    /* Destruímos el plan */  
    fftw_destroy_plan(p);  
 
    /* Liberamos memoria */  
    free(out);  
    free(in);  
  }  
  return 0;  
}  
 
int compute_number_of_samples(FILE *input_file) {  
  fseek(input_file,0,SEEK_END);  
  int samples = ftell(input_file)/sizeof(float);  
  rewind(input_file);  
  return samples;  
}  
 
read_signal(double *in, FILE *input_file, int samples) {  
  int i;  
  for(i=0; i<samples; i++) {  
    float x;  
    fread(&x,sizeof(float),1,input_file);  
    in[i] = (double)x;  
    spin();  
  }  
}