/*  
 * demodulator.c -- Desmodula una señal (desplaza su espectro).  
 *  
 * Este fichero fuente puede encontrarse en:  
 * http://www.ace.ual.es/~vruiz/docencia/redes/practicas/demodulator.c  
 *  
 * Compilar escribiendo (el paquete fftw debería estar instalado!):  
 * gcc demodulator.c -o demodulator 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 carrier_frequency modulated_signal.float > unmodulated_signal.float\n"  
            ,argv[0]);  
  } else {  
 
    FILE *input_file = fopen(argv[2],"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       *signal   = (double       *)fftw_malloc(samples*sizeof(double      ));  
    fftw_complex *spectrum = (fftw_complex *)fftw_malloc(samples*sizeof(fftw_complex));  
    fftw_complex *tmp      = (fftw_complex *)fftw_malloc(samples*sizeof(fftw_complex));  
 
    /* Creamos un plan para la transformada directa en inversa */  
    fftw_plan f = fftw_plan_dft_r2c_1d(samples, signal, spectrum, 0);  
    fftw_plan b = fftw_plan_dft_c2r_1d(samples, spectrum, signal, 0);  
 
    read_signal(signal, input_file, samples);  
 
    /* Calculamos la transformada de Fourier */  
    fftw_execute(f);  
 
    restore_spectrum(spectrum, tmp, samples, atoi(argv[1]));  
 
    /* Calculamos la transformada de Fourier inversa */  
    fftw_execute(b);  
 
    write_signal(signal, stdout, samples);  
 
    /* Destruímos los planes */  
    fftw_destroy_plan(f); fftw_destroy_plan(b);  
 
    /* Liberamos memoria */  
    free(tmp);  
    free(spectrum);  
    free(signal);  
  }  
  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 *signal, FILE *input_file, int samples) {  
  int i;  
  for(i=0; i<samples; i++) {  
    float x;  
    fread(&x,sizeof(float),1,input_file);  
    signal[i] = (double)x;  
    spin();  
  }  
}  
 
write_signal(double *signal, FILE *output_file, int samples) {  
  int i;  
  for(i=0; i<samples; i++) {  
    float sample = signal[i]/samples;  
    fwrite(&sample,sizeof(float),1,output_file);  
    spin();  
  }  
}  
 
restore_spectrum(fftw_complex *spectrum, fftw_complex *tmp, int samples, int carrier) {  
  int i;  
  for(i=0; i<samples/2; i++) {  
    tmp[i][0] = 0;  
    tmp[i][1] = 0;  
    spin();  
  }  
  for(i=0; i<samples/2 - carrier; i++) {  
    tmp[i][0] = spectrum[carrier+i][0];  
    tmp[i][1] = spectrum[carrier+i][1];  
    spin();  
  }  
  for(i=0; i<samples/2; i++) {  
    spectrum[i][0] = tmp[i][0];  
    spectrum[i][1] = tmp[i][1];  
    spin();  
  }  
}