#include <iostream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <complex>
#include "AudioFile.h"

/*
  xmgrace Transformed.dat

*/

using namespace std;

int main()
{
    //---------------------------------------------------------------
    // 1. Set a file path to an audio file on your machine

    const std::string base   = "sine-wave";
    const std::string input  = base+".wav";

    //---------------------------------------------------------------
    // 2. Create an AudioFile object and load the audio file

    AudioFile<float> a;
    bool loadedOK = a.load(input);
    assert (loadedOK);

    cout << "Read: "              << input                  << endl;
    cout << "Bit Depth: "         << a.getBitDepth()        << endl;
    cout << "Sample Rate: "       << a.getSampleRate()      << endl;
    cout << "Num Channels: "      << a.getNumChannels()     << endl;
    cout << "Length in Seconds: " << a.getLengthInSeconds() << endl;

    double rate=a.getSampleRate();
    int blocksize = 4410;
    int nblocks = a.getNumSamplesPerChannel()/blocksize;

    double fmin=rate/blocksize;
    double fmax=rate/2;

    cout << "Number of blocks: "                   << nblocks               << endl;
    cout << "Period of one block: "                << blocksize/rate << "s" << endl;

    cout << "Minimal frequency due to block size " << fmin << " Hz" << endl;
    cout << "Maximal frequency due sampling      " << fmax << " Hz" << endl;
    //---------------------------------------------------------------
    // 3. Let's Fourier transform all the blocks

    int channel=0;
    vector<complex<float> > amplitudesstore(nblocks*blocksize, 0.0);
    vector<complex<float> > e(blocksize*blocksize);

    //---------------------------------------------------------------
    // 4. Lets precalculate exp(i 2pi j*k / blocksize)   for speed

    for (int k=0; k < blocksize ; k++)
        for (int j=0 ; j < blocksize; j++)
        {
            int o=k*blocksize+j;
            complex<double> z(0.0, 2.0*M_PI*j*k/blocksize);
            e[o] = exp(z);
        }

    //---------------------------------------------------------------
    // 5. Lets calculate fourier transform the individual blocks

    for (int b = 0; b < nblocks; b++)
        for (int k=0; k < blocksize ; k++)
        {
            int offset=b*blocksize;
            complex<double> ak(0.0,0.0);

            // a_k = sum_j=0 ^blocksize-1  f(j) exp(I 2pi j k j/N)
            for (int j=0 ; j < blocksize; j++)
                ak += e[k*blocksize+j] * a.samples[channel][offset+j];

            // Store ak
            amplitudesstore[offset+k]=ak;
        }


    //
    // Here we could transform the amplitudes in some interesting way.

     for (int b = 0; b < nblocks; b++)
        for (int k=0; k < blocksize ; k++)
           {
               double f=0.0;
               if (k<blocksize/2) f=k*fmin;              // positive
                             else f=(k-blocksize)*fmin;  // Negative


               // Kill all frequencies below 300hz, noting both positive and
               // negative frequencies.

               if (fabs(f)<300) amplitudesstore[b*blocksize+k]=complex(0.0,0.0);
           }

    //---------------------------------------------------------------
    // 6. Lets calculate exp(-i 2pi j*k / blocksize)

    for (int k=0; k < blocksize ; k++)
        for (int j=0 ; j < blocksize; j++)
        {
            int o=k*blocksize+j;
            complex<double> z(0.0, -2.0*M_PI*j*k/blocksize);
            e[o] = exp(z);
        }

    //---------------------------------------------------------------
    // 7. Lets calculate the inverse fourier transform the individual blocks

    for (int b = 0; b < nblocks; b++)
        for (int j=0; j < blocksize ; j++)
        {
            complex<double> signal(0.0,0.0);

            // f(j) = sum_k=0 ^blocksize-1  a_k exp(-I 2pi k j/N)
            for (int k=0 ; k < blocksize; k++)
                signal += amplitudesstore[b*blocksize+k]*e[k*blocksize+j];

            a.samples[channel][b*blocksize+j]=signal.real()/blocksize;
        }

    const std::string output  = base+"_transformed.wav";
    const std::string output2 = base+"_transformed.dat";

    a.save (output, AudioFileFormat::Aiff);

    ofstream fo(output2);
    for (int i = 0; i < a.getNumSamplesPerChannel(); i++)
    {
        double t=1.0*i/rate;
        fo <<  t << " " << a.samples[channel][i] << endl;
    }
    fo.close();
}