#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
#define N 1764                  //size of array input to FFT

int main()
{
    double *data = (double *) malloc(sizeof(double) * N);       //input data of size n
    fftw_complex *fft = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * N); //fftw data of size n
    double x, cc, ts;           //x is a dummy variable, cc is the output*complex conjugate of the FFT output, ts is the time of each sample
    fftw_plan plan;
    FILE *voice;
    FILE *out;
    FILE *ccout;
    FILE *time;
    int i, j, p;
    int ch;
    long nxl = 0;
    FILE *fp;
    double a, mx;

    ts = N / 44100;             //number of samples taken/number of samples per second

    if ((fp = fopen("sinwav.dat", "r")) == NULL) {        //open the file for reading
        printf("Can`t open sinwav.dat\n");
        exit(1);
    }


    while ((ch = getc(fp)) != EOF) {     //count the number of line
        if (ch == '\n')
            nxl++;
    }

    fclose(fp);                 //close this

    printf("File sinwav.dat has %ld lines\n", nxl);     //the number of lines in the file
    time = fopen("time.dat", "w");      //create the file for writing to
    if ((voice = fopen("anna.dat", "r")) == NULL) {
        printf("Cannot open %s\n", "anna.dat");
        exit(EXIT_FAILURE);
    }

    for (p = 0; p < nxl / N; p++) {       //repeat the following for the number of lines in file divided by the number of files per input to FFT

        for (i = 0; i < N; i++) {       //input lines from file into data
            fscanf(voice, "%lf %lf", &x, &data[i]);     //I think I'm right in saying fscanf carries on scanning from the file stream, so will pick up where it left off
        }

        plan = fftw_plan_dft_r2c_1d(N, data, fft, FFTW_ESTIMATE);       //plan the FFT
        fftw_execute(plan);     //execute it

        for (j = 0; j < N; j++) {       //find the maximum value for a given FFT

            cc = fft[0][j] * fft[0][j] + fft[1][j] * fft[1][j];

            if (j = 0) {
                mx = (j * 44100) / (2 * N);
                a = cc;
            } else {

                if (cc > a) {
                    mx = (j * 44100) / (2 * N);
                    a = cc;
                }
            }
        }
        fprintf(time, "%lf\t%lf\t%lf\n", p * ts, mx, a);        //print this to file
    }
    fclose(out);
    fclose(ccout);
    fclose(time);
    exit(0);
}