// これはメインDLL ファイルです。

#include "myFFT.h"

#include <math.h>

//*************************************************************************
#define _MY_MAKE_DLL_
//*************************************************************************

/*
Fast Fourier/Cosine/Sine Transform
    dimension   :one
    data length :power of 2
    decimation  :frequency
    radix       :split-radix
    data        :inplace
    table       :use
functions
    cdft: Complex Discrete Fourier Transform
    rdft: Real Discrete Fourier Transform
    ddct: Discrete Cosine Transform
    ddst: Discrete Sine Transform
    dfct: Cosine Transform of RDFT (Real Symmetric DFT)
    dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
function prototypes
    void cdft(int, int, double *, int *, double *);
    void rdft(int, int, double *, int *, double *);
    void ddct(int, int, double *, int *, double *);
    void ddst(int, int, double *, int *, double *);
    void dfct(int, double *, double *, int *, double *);
    void dfst(int, double *, double *, int *, double *);
macro definitions
    USE_CDFT_PTHREADS : default=not defined
        CDFT_THREADS_BEGIN_N  : must be >= 512, default=8192
        CDFT_4THREADS_BEGIN_N : must be >= 512, default=65536
    USE_CDFT_WINTHREADS : default=not defined
        CDFT_THREADS_BEGIN_N  : must be >= 512, default=32768
        CDFT_4THREADS_BEGIN_N : must be >= 512, default=524288

-------- Complex DFT (Discrete Fourier Transform) --------
    [definition]
        <case1>
            X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
        <case2>
            X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
        (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
    [usage]
        <case1>
            ip[0] = 0; // first time only
            cdft(2*n, 1, a, ip, w);
        <case2>
            ip[0] = 0; // first time only
            cdft(2*n, -1, a, ip, w);
    [parameters]
        2*n            :data length (int)
                        n >= 1, n = power of 2
        a[0...2*n-1]   :input/output data (double *)
                        input data
                            a[2*j] = Re(x[j]), 
                            a[2*j+1] = Im(x[j]), 0<=j<n
                        output data
                            a[2*k] = Re(X[k]), 
                            a[2*k+1] = Im(X[k]), 0<=k<n
        ip[0...*]      :work area for bit reversal (int *)
                        length of ip >= 2+sqrt(n)
                        strictly, 
                        length of ip >= 
                            2+(1<<(int)(log(n+0.5)/log(2))/2).
                        ip[0],ip[1] are pointers of the cos/sin table.
        w[0...n/2-1]   :cos/sin table (double *)
                        w[],ip[] are initialized if ip[0] == 0.
    [remark]
        Inverse of 
            cdft(2*n, -1, a, ip, w);
        is 
            cdft(2*n, 1, a, ip, w);
            for (j = 0; j <= 2 * n - 1; j++) {
                a[j] *= 1.0 / n;
            }
 //以下略、大浦様のプログラムをコピーして下さい。       .

//************************************************************************************************************

//ここから、
//************************************************************************************************************

   double pi = 3.14159265358979;

double func_Rect(int i,int n)
{
	return 1.0;//Rectangular
}

double func_Ham(int i,int n)
{
	//return 0.54-0.46*cos(2*pi*i/(N-1));//Hamming
	return (double)( 0.54-0.46*cos(2*pi*i/n) );//Hamming

}

double func_Han(int i,int n)
{
	//return 0.5-0.5*cos(2*pi*i/(N-1));//Hanning
	return (double)( 0.5-0.5*cos(2*pi*i/n) );//Hanning

}

double func_Blkman(int i,int n)
{
	//return 0.42-0.5*cos(2*pi*i/(N-1))+0.08*cos(4*pi*i/(N-1));//Blackman
	return (double)( 0.42-0.5*cos(2*pi*i/(n-1))+0.08*cos(4*pi*i/n) );//Blackman

}

double func_BlkmanHarris(int i,int n)//Blackman-Harris
{
	//return 0.35875-0.48829*cos(2*pi*i/(N-1))+0.14128*cos(4*pi*i/(N-1))-0.01168*cos(6*pi*i/(N-1));
	return (double)( 0.35875-0.48829*cos(2*pi*i/n)+0.14128f*cos(4*pi*i/n)-0.01168f*cos(6*pi*i/n) );

}

double func_sin(int i,int n)//sin
{
	return (double)( sin(pi*(i+0.5)/(2.0*n)) );
}

typedef double (*mfunc)(int,int);

mfunc Func[]={func_Rect,func_Han,func_Ham,func_Blkman,func_BlkmanHarris,func_sin};

void WindowFunc(int func_no,int number_of_spectrum,double* source,double* result)
{
	mfunc myFunc=Func[func_no];

	for (int i=0;i<number_of_spectrum;i++)
	{
		
		*result++= (double)(*source++ * myFunc(i,number_of_spectrum));
	}
	return ;
}
void inv_WindowFunc(int func_no,int number_of_spectrum,double* source,double* result)
{
	mfunc myFunc=Func[func_no];

	for (int i=0;i<number_of_spectrum;i++)
	{
		
		*result++= (double)(*source++ / myFunc(i,number_of_spectrum));
	}
	return ;
}
//**************************************************************************************************
//Fir fliter class
//**************************************************************************************************

        //窓関数法による、fir filterの作成　

        //returnは、filter_size
double* CFir:: get_LPFir_coeff( int selected_window, int sample_per_sec, int size, int edge_f )
        {
			//フィルターサイズを指定して、遷移周波数は、計算で求める
		fc = 0;     //cut off frequency
        fe = 0;     //edge frequency
        
		filter_size = size;
        fe = edge_f;

            // use window
            switch (selected_window)
            {
                case 0: delta_frequency_current = (double)(0.9 / filter_size); break;
                case 1: delta_frequency_current = (double)(3.1 / filter_size); break;
                case 2: delta_frequency_current = (double)(3.3 / filter_size); break;
                case 3: delta_frequency_current = (double)(5.5 / filter_size); break;
                case 4:
                case 5:
                default: delta_frequency_current = (double)(3.3 / filter_size); break;
            }

            delta_frequency_current *= sample_per_sec;//de_normalize

            fc = fe + delta_frequency_current / 2.0;    // not normalize

            //normalize
            fc /= sample_per_sec;
//don't foreget to delete
			return_array = new double[filter_size];

            int k = 0;
            for (int i = -filter_size / 2; i < filter_size / 2; i++)//LPFirの係数をずらす。
            {
                if (i == 0)
                    return_array[k] = 2.0 * fc;
                else
                {
                    return_array[k] = 2.0 * fc * sin(2 * pi * fc * i) / (2 * pi * fc * i);
                }

                //taps_current[K] *= Math.Pow(2, 6);  //scaling 規定値
                //return_array[K] *= Math.Pow(2, 2);  //scaling

                k++;
                //if (K == filter_size / 2) break;
            }

            return return_array;

        }
        // double LPFir_coef_FFT_result;

double* CFir::get_windowed_LPFir_coef_FFT_result( double* coef, int coef_length,int capture_data_length,int func_no, 
																			int number_of_spectrum )
{
	int i;
//don't foreget to delete
	return_array2 = new double[number_of_spectrum];

	if( coef_length +capture_data_length -1 >= number_of_spectrum )
		return (double*)0;

	//ここで、直線たたみ込みのために、ゼロ詰めする。filter coef length ->numberOfspectrum
	// x(N),h(M) とすると　numberOfspectrum >= N+M-1
    for( i = 0; i < coef_length; i++ )
		return_array2[i]=coef[i];

	for( ;i < number_of_spectrum;i++ )
		return_array2[i]=0;

	//ここで、窓関数を掛ける
	WindowFunc( func_no, number_of_spectrum, return_array2 ,return_array2);

	//cdft(2*numberOfspectrum,-1,2*retrun_array2,ip,w);
	coef_fft  = new double[ 2 * number_of_spectrum ];
	
	is_first_time=true;

	if(is_first_time)
	{
		ip = new int[ (int)(2 + sqrt( (double)number_of_spectrum )) ];
		ip[0]=0;
		w  = new double[ number_of_spectrum / 2 ];
		is_first_time=false;
	}

	for(i=0;i<number_of_spectrum;i++)
	{
		coef_fft[2*i]=return_array2[i];
		coef_fft[2*i+1]=0;
	}
	// do FFT
	cdft(2*number_of_spectrum,-1,coef_fft,ip,w);

	return coef_fft;

}

double* CFir::get_windowed_CaptureData_FFT_result(float* capture_data,int capture_data_length, int func_no,
												int coef_length, int number_of_spectrum )
{
	int i;
	//フィルター係数coefの次数は、入力配列の次数より、小さい事
	if(capture_data_length + coef_length -1 >=number_of_spectrum)
		return (double*)0;

	if(is_first_time)
		return_array2=new double[number_of_spectrum];

	//ゼロ詰め。巡回たたみこみを避けるために

	for(i=0;i<capture_data_length;i++)
		return_array2[i]=(double)capture_data[i];
	for(;i<number_of_spectrum;i++)
		return_array2[i]=0.0;

	//ここで、窓関数を掛ける
	WindowFunc( func_no, number_of_spectrum, return_array2 ,return_array2);
	
	//FFTの準備
	capture_fft=new double[2*number_of_spectrum];
	
	for(i=0;i<number_of_spectrum;i++)
	{
		capture_fft[2*i]=return_array2[i];
		capture_fft[2*i+1]=0.0;
	}

	if(is_first_time)
	{
		ip=new int[(int)(2+sqrt((double)number_of_spectrum))];
		ip[0]=0;
		w=new double[number_of_spectrum/2];
		is_first_time=false;
	}
	//do FFT
	cdft(2*number_of_spectrum,-1,capture_fft,ip,w);

	return capture_fft;
	
}

//input ='src'= src1,'dest'= src2 : out= 'dest' pos
//return length = 2*number_of_spectrum
double*	CFir::multi_FFT_results( double* src, double* dest, int number_of_spectrum )
{
	int i;
	/*for(i=0;i<2*number_of_spectrum;i++)//これは間違い	(a+bi)*(c+di)=ac-bd+(ad+bc)i
		dest[i]=src[i]*dest[i];*/

	//正しくは
	for(i=0;i<number_of_spectrum;i++)
	{
      	 dest[2*i]   = src[2*i] * dest[2*i]   - src[2*i+1] * dest[2*i+1];
             dest[2*i+1] = src[2*i] * dest[2*i+1] + src[2*i+1] * dest[2*i];	
	}
	
	return dest;
}

float* CFir::get_data_from_mulit_FFT_results( double* src, int number_of_spectrum, int capturedata_length, 
										   int func_no )
{
	int i;
	double* result=new double[number_of_spectrum];

	//IFFTの準備
	if(is_first_time)
	{
		ip = new int[ (int)(2 + sqrt( (double)number_of_spectrum )) ];
		ip[0]=0;

		w  = new double[ number_of_spectrum / 2 ];
		is_first_time=false;
	}	
	//do IFFT
	cdft(2*number_of_spectrum,1,src,ip,w);
	
	for(i=0;i<2*number_of_spectrum-1;i++)
	{
		/src[i] *=1.0/number_of_spectrum;

	}

	//実部のみ、取り出す
	for(i=0;i<number_of_spectrum;i++)
		result[i]=src[2*i];

	//逆窓掛け
	inv_WindowFunc( func_no, number_of_spectrum, result , result);

	//CaptureData_lengthだけ、実部を取り出す
	returned_data = new float[capturedata_length];
	for(i=0;i<capturedata_length;i++)
		returned_data[i] = (float)result[i];
	
	delete[] result;
	return returned_data;

}

CFir::~CFir()
{
	//if( return_array2 != (void*)0 )
		delete[] return_array2;
	//if( return_array != (void*)0 )
		delete[] return_array;
	//if( coef_fft != (void*)0)
	//	delete[] coef_fft;
	//if( capture_fft !=(void*)0)
	//	delete[] capture_fft;
	//if( ip != (void*)0)
		delete[] ip;
	//if( w != (void*)0)
		delete[] w;
	//if(returned_data !=(void*)0)
		delete[] returned_data; 
}
CFir::CFir()
{
}