#ifndef WIN4C_H__
#define WIN4C_H__
#pragma once

#define PI  3.141592653589793

#define EPS  1e-8

#include <vector>

std::vector<double> rectangle(int);
std::vector<double> boxcar(int); //same as matlab, rectangle window.
std::vector<double> bartlett(int);
std::vector<double> hanning(int, bool symmetric=true);
std::vector<double> hamming(int);
std::vector<double> blackman(int);
std::vector<double> kaiser(int, double alpha=0.5);
std::vector<double> gausswin(int, double alpha=2.5);
std::vector<double> chebwin(int, double r = 100);
std::vector<double> triang(int);
std::vector<double> tukeywin(int, double r = 0.5);
double I0(double alpha);

/*
   NOTE: Use the HANN function to get a Hanning window which has the
          first and last zero-weighted samples.
*/
std::vector<double> hann(int, bool symmetric = true);

#endif

#include "win4c.h"
#include "fftw3.h"


std::vector<double> boxcar(int n)
{
    return rectangle(n);
}

std::vector<double> rectangle(int n)
{
    std::vector<double> win(n, 1);

    return win;
}

std::vector<double> bartlett(int n)
{
    std::vector<double> win(n);
    for (int i = 0; i < (n + 1) / 2; ++i)
    {
        win[i] = 2.0 * i / (n - 1.0);
        win[n - 1 - i] = win[i];
    }
    return win;
}

std::vector<double> hanning(int n, bool symmetric)
{
    std::vector<double> win(n);
    if (symmetric)
    {
        for (int i = 0; i < (n + 1) / 2; ++i)
        {
            win[i] = 0.5 - 0.5 * cos(2 * PI * (i + 1) / (n + 1));
            win[n - 1 - i] = win[i];
        }
    }
    else
    {
        win.push_back(0);
        std::vector<double> tmp = hanning(n - 1, true);
        win.insert(win.begin() + 1, tmp.begin(), tmp.end());
    }
    return win;
}

std::vector<double> hamming(int n)
{
    std::vector<double> win(n);
    for (int i = 0; i < (n + 1) / 2; ++i)
    {
        win[i] = (0.54 - 0.46 * cos (2 * PI * i / (n - 1.0)));
        win[n - 1 - i] = win[i];
    }
    return win;
}

std::vector<double> blackman(int n)
{
    std::vector<double> win(n);
    for (int i = 0; i < (n + 1) / 2; ++i)
    {
        win[i] = (0.42 - 0.50 * cos(2 * PI * i / (n - 1.0)) + 
                  0.08 * cos(2 * 2 * PI * i / (n - 1.0)));
        win[n - 1 - i] = win[i];
    }
    return win;
}

std::vector<double> kaiser(int n, double alpha)
{
    std::vector<double> win(n);
    for (int i = 0; i < (n + 1) / 2; ++i)
    {
        double beta = 2 * alpha * sqrt(i*(n - i - 1.0)) / (n - 1.0);
        win[i] = I0(beta) / I0(alpha);
        win[n - 1 - i] = win[i];
    }
    return win;
}

std::vector<double> gausswin(int n, double alpha)
{
    std::vector<double> win(n);
    double center = (n - 1) / 2.0;

    for (int i = 0; i < (n + 1) / 2; ++i)
    {
        double tmp = alpha * (i - center) / center;
        win[i] = exp(-0.5 * tmp * tmp);
        win[n - 1 - i] = win[i];
    }
    return win;
}

std::vector<double> chebwin(int M, double r)
{
    std::vector<double> win(M);
    double order = M - 1.0;
    double beta = cosh(1.0 / order * acosh(pow(10, (abs(r) / 20.))));
    double *x = new double[M];
    for (int i = 0; i < M; i++)
    {
        double temp = cos(PI * i / M)*beta;
        if (temp > 1)
        {
            temp = cosh(order * acosh(temp));
        }            
        else if (temp < -1)
        {
            temp = (2 * (M % 2) - 1) * cosh(order * acosh(-temp));
        }
        else
        {
            temp = cos(order * acos(temp));
        }            
        x[i] = temp;
    }
    fftw_complex *fftIn = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * M);
    fftw_complex *fftOut = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * M);
    fftw_plan plan = fftw_plan_dft_1d(M, fftIn, fftOut, FFTW_FORWARD, FFTW_MEASURE);
    double maxvalue(-1e8);

    if (M % 2)//M ����
    {
        for (int j = 0; j < M; j++)
        {
            fftIn[j][0] = x[j];
            fftIn[j][1] = 0;
        }
        fftw_execute(plan);
        int n = (M + 1) / 2; 
        win[n - 1] = fftOut[0][0];
        maxvalue = fftOut[0][0];
        for (int j = n - 1; j >= 1; j--)
        {
            win[n - 1 - j] = fftOut[j][0];
            win[j + n - 1] = fftOut[j][0];
            if (maxvalue < fftOut[j][0])
            {
                maxvalue = fftOut[j][0];
            }                
        }
    }
    else
    {
        for (int j = 0; j < M; j++)
        {
            fftIn[j][0] = x[j] * cos(PI / M * j);
            fftIn[j][1] = x[j] * sin(PI / M * j);
        }
        fftw_execute(plan);
        int n = M / 2 + 1;
        for (int j = n - 1; j >= 1; j--)
        {
            win[n - 1 - j] = fftOut[j][0];
            win[j + n - 2] = fftOut[j][0];
            if (maxvalue < fftOut[j][0]) 
            {
                maxvalue = fftOut[j][0];
            }                
        }
    }
    for (int j = 0; j < M; j++)
    {
        win[j] = win[j] / maxvalue;
    }
    delete [] x;
    x = nullptr;
    fftw_destroy_plan(plan);
    fftw_free(fftIn);
    fftw_free(fftOut);
    return win;
}

std::vector<double> triang(int n)
{
    std::vector<double> win(n);
    for (int i = 0; i < (n + 1) / 2; ++i)
    {
        win[i] = (2.0 * ( i + 1) - 1 + n % 2) / (n + n % 2);
        win[n - 1 - i] = win[i];
    }
    return win;
}

std::vector<double> hann(int n, bool symmetric)
{
    std::vector<double> win(n);
    win.push_back(0);

    std::vector<double> tmp = hanning(n - 2, symmetric);
    win.insert(win.begin() + 1, tmp.begin(), tmp.end());
    win.push_back(0);
    return win;
}

std::vector<double> tukeywin(int n, double r)
{
    if (r <= 0)
    {
        return boxcar(n);
    }
    if (r >= 1)
    {
        return hann(n);
    }
    std::vector<double> win(n);

    std::vector<double> t(n);
    double step = 1.0 / (n - 1);
    for (int i = 0; i < n - 1; i++)
    {
        t[i] = i * step;
    }
    t[n - 1] = 1;

    double per = r / 2.0;
    int tl = floor(per * (n - 1)) + 1;
    int th = n - tl + 1;
    //Window is defined in three sections : taper, constant, taper
    //    w = [((1 + cos(pi / per * (t(1:tl) - per))) / 2);  ones(th - tl - 1, 1); ((1 + cos(pi / per * (t(th:end) - 1 + per))) / 2)];
    for (int i = 0; i < tl; i++)
    {
        win[i] = (1 + cos(PI / per * (t[i] - per))) / 2.0;
    }
    for (int i = tl; i < th - 1; i++)
    {
        win[i] = 1;
    }
    for (int i = th - 1; i < n; i++)
    {
        win[i] = (1 + cos(PI / per * (t[i] - 1 + per))) / 2.0;
    }
    return win;
}
double I0(double alpha)
{
    double dNew;
    double K = alpha / 2.0;
    const int MAXTERM = 25 + 1;
    double J = 1.0, dOld = 1.0;
    bool converge = false;

    // Use series expansion definition of Bessel.
    for (int i = 1; i < MAXTERM; ++i)
    {
        J *= K / i;
        dNew = dOld + J * J;

        if ((dNew - dOld) < EPS)
        {
            converge = true;
            break;
        }
        dOld = dNew;
    }

    if (!converge)
        return 0;

    return dNew;
}