/* ft computer.
   Usage:
   ndft freq_beg freq_end delta_freq < input > output

   reads input data from standard input
     first line of input is ignored
     (normally this line contains the number of points in file)
   writes output on standard output
*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define MAX_LINE 	80
#define MAX_NAME	120
#define TWOPI		6.283185307
#define	NPTOS		20000
int main(int argc,char *argv[])
{
double	t[NPTOS]	,
	y[NPTOS]	;
double	t_inicio	,
	avg		;
double	f_beg	,
	f_end	,
	f_step	,
	f	;
int	i	,
	npts	;
FILE	*fp_inp	,
	*fp_out	;
char	data_file[MAX_NAME]	,
	sp_file[MAX_NAME]	,
	lixo[MAX_LINE]		;


f_beg = atof(argv[1])	;
f_end = atof(argv[2])	;
f_step= atof(argv[3])	;

/* reads the data file skipping first line */
 for (i=0 ; (scanf("%lf%lf",&t[i],&y[i]) != EOF) ; i++)	;
 npts = i	;

/* makes initial time = 0, computes average */
t_inicio = t[0]	;
for (i=0 ; i<npts ; i++)
	{
	avg  = avg  + y[i]	;
	t[i] = t[i] - t_inicio	;
	}
avg = avg / npts	;

/* average = 0 */

for (i=0 ; i<npts ; i++)
	y[i] = y[i] - avg ;

/* starts computing FT */

/* most external loop in frequency */
 for (f=f_beg ; f<=f_end ; f += f_step)
   {
     double	fr = 0.	,
       fi = 0.	,
       amplitude=0.0	;
     double	sqrt()	;
     /* internal loop in data points */
     for (i=0 ; i<npts ; i++)
       {
	 double 	arg	,
	   c	,
	   s	;
	 arg = TWOPI * f * t[i]	;
	 c   = cos(arg)		;
	 s   = sin(arg)		;
	 fr += y[i] * c		;
	 fi += y[i] * s		;
         // printf("%lf %lf %lf %lf %lf %lf %lf\n",t[i],f,arg,c,fr,s,fi);
       }
     fr /= npts	;
     fi /= npts 	;
     amplitude = fr*fr + fi*fi	;
     amplitude = sqrt(amplitude)	;
     printf("%le\t%le\n",f,amplitude)	;
   }
}
