using System;


namespace WindowsFormsApplication1
{
    public class Interpolator
    {


        //-------------------------------------------------------------------------------------------------
        //-------------------------DFT for counting the spectrum of non power of two array-----------------
        //-------------------------------------------------------------------------------------------------
        private static void dft(double[] inreal, double[] inimag, double[] outreal, double[] outimag)
        {
            int n = inreal.Length;
            for (int k = 0; k < n; k++)
            {  // For each output element
                double sumreal = 0;
                double sumimag = 0;
                for (int t = 0; t < n; t++)
                {  // For each input element
                    sumreal += inreal[t] * Math.Cos(2 * Math.PI * t * k / n) + inimag[t] * Math.Sin(2 * Math.PI * t * k / n);
                    sumimag += -inreal[t] * Math.Sin(2 * Math.PI * t * k / n) + inimag[t] * Math.Cos(2 * Math.PI * t * k / n);
                }
                outreal[k] = sumreal;
                outimag[k] = sumimag;
            }
        }
        //-------------------------------------------------------------------------------------------------
        //  Standard Cooley-Tukey FFT algorithm, based on routine of N. Brenpairser, as described
        //  in "Numerical Recipes"
        //
        //  Note that the data file for this routine is standard 0 based C type routine.  NOT the
        //  FORTRAN style 1 based indexing used in "Numerical recipes in C"

        /// <summary>
        /// Cooley-Tukey fast Fourier transform, zero index (C style) based
        /// </summary>
        /// <param name="data">the data are alternating real/imaginary pairs</param>
        /// <param name="npairs">the number of real/imaginary pairs (size of data[]/2)</param>
        /// <param name="sign">1 for forward transform, -1 for inverse transform</param>
        private static void fft(double[] data, int npairs, int sign)
        {
            int n, mmax, m, j, istep, i;
            double wtemp, wr, wpr, wpi, wi, theta;
            double temp, tempr, tempi;

            const double two_Pi = 2.0 * Math.PI;

            // put data in bit-reversed order 
            n = 2 * npairs;
            j = 1;

            for (i = 1; i < n; i += 2)
            {
                if (j > i)
                {
                    temp = data[j - 1]; data[j - 1] = data[i - 1]; data[i - 1] = temp;
                    temp = data[j]; data[j] = data[i]; data[i] = temp;
                }

                m = n / 2;

                while (m >= 2 && j > m)
                {
                    j = j - m;
                    m = m / 2;
                }

                j = j + m;
            }

            mmax = 2;

            while (n > mmax)
            {
                istep = 2 * mmax;
                theta = two_Pi / (double)(-sign * mmax);
                wtemp = Math.Sin(0.5 * theta);
                wpr = -2.0 * wtemp * wtemp;
                wpi = Math.Sin(theta);
                wr = 1.0;
                wi = 0.0;
                for (m = 1; m < mmax; m += 2)
                {
                    for (i = m; i <= n; i += istep)
                    {
                        j = i + mmax;
                        tempr = wr * data[j - 1] - wi * data[j];
                        tempi = wr * data[j] + wi * data[j - 1];
                        data[j - 1] = data[i - 1] - tempr;
                        data[j] = data[i] - tempi;
                        data[i - 1] += tempr;
                        data[i] += tempi;
                    }
                    wtemp = wr;
                    wr = wtemp * wpr - wi * wpi + wr;
                    wi = wi * wpr + wtemp * wpi + wi;
                }
                mmax = istep;
            }

            /*  if (sign == 1) 
              {
                  for (n = 0; n < npairs; n+=2)
                  {
                      data[n] = data[n] / npairs;
                      data[n+1] = data[n+1] / npairs;
                  }
              }*/
        }

        //-------------------------------------------------------------------------------------------------
        //-------------------------                Interpolation    ---------------------------------------
        //-------------------------------------------------------------------------------------------------
        public static void Interpolate(int n_times, byte[] origdata, byte[] intpdata)
        {
            int n;

            int len;
            double comp;
            double from;

            double[] DFT_inreal;
            double[] DFT_inimag;
            double[] DFT_outreal;
            double[] DFT_outimag;

            double[] intpdata_big = new double[1024 * 2];// *2 becouse they will be complex



            switch (n_times)
            {
                case 2: //2x
                    // in final 533 points are needed. Nearest bigger 2^n is 1024
                    // so 1024/2 = 512
                    len = 512;
                    comp = 512;
                    from = 266.0 - len / 2.0;

                    double[] intpdata_spect = new double[512 * 2];// *2 becouse they will be complex

                    // taking appropirate samples from mainbuff, 
                    //preparing to create spectrum
                    for (n = 0; n < len; n++)
                    {
                        intpdata_spect[n * 2] = origdata[(int)(n + from)];
                    }

                    // FFT 
                    fft(intpdata_spect, 512, 1);




                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 0; n < 512; n++)
                    {
                        intpdata_big[n] = intpdata_spect[n] / comp; ;

                    }

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 513; n < intpdata_spect.Length; n++)
                    {
                        intpdata_big[n + 1024] = intpdata_spect[n] / comp; ;

                    }

                    // IFFT of widened spectrum
                    fft(intpdata_big, 1024, -1);

                    for (n = 0; n < 533; n++) //copy convert and saturate
                    {
                        if (intpdata_big[(n + 512 - 266) * 2] > 255.0)
                            intpdata[n] = 255;
                        else if (intpdata_big[(n + 512 - 266) * 2] < 0.0)
                            intpdata[n] = 0;
                        else
                            intpdata[n] = (byte)intpdata_big[(n + 512 - 266) * 2];
                                
                    }

                    break;


                case 5:   //5x

                    // in final 533 points are needed. Nearest bigger 2^n is 1024
                    // so 1024/5 = 204.8
                    len = 205;
                    comp = 204.8;
                    from = 266.0 - len / 2.0;

                    DFT_inreal = new double[len];
                    DFT_inimag = new double[len];
                    DFT_outreal = new double[len];
                    DFT_outimag = new double[len];


                    // taking appropirate samples from mainbuff, 
                    //preparing to create spectrum
                    for (n = 0; n < len; n++)
                    {
                        DFT_inreal[n] = origdata[(int)(n + from)];
                    }

                    // non 2^n spectrum
                    dft(DFT_inreal, DFT_inimag, DFT_outreal, DFT_outimag);

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 0; n < 101 * 2 + 1; n += 2)
                    {
                        intpdata_big[n] = DFT_outreal[n / 2] / comp;
                        intpdata_big[n + 1] = DFT_outimag[n / 2] / comp;
                    }

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 104 * 2; n < DFT_outreal.Length * 2; n += 2)
                    {
                        intpdata_big[n + 819 * 2] = DFT_outreal[n / 2] / comp;
                        intpdata_big[n + 1 + 819 * 2] = DFT_outimag[n / 2] / comp;
                    }

                    // IFFT of widened spectrum
                    fft(intpdata_big, 1024, -1);

                    for (n = 0; n < 533; n++) //copy convert and saturate
                    {
                        if (intpdata_big[(n + 512 - 266) * 2] > 255.0)
                            intpdata[n] = 255;
                        else if (intpdata_big[(n + 512 - 266) * 2] < 0.0)
                            intpdata[n] = 0;
                        else
                            intpdata[n] = (byte)intpdata_big[(n + 512 - 266) * 2];

                    }

                    break;


                case 10:   //10x

                    // in final 533 points are needed. Nearest bigger 2^n is 1024
                    // so 1024/10 = 102.4
                    len = 102;
                    comp = 102.4;
                    from = 266.0 - len / 2.0;

                    DFT_inreal = new double[len];
                    DFT_inimag = new double[len];
                    DFT_outreal = new double[len];
                    DFT_outimag = new double[len];


                    // taking appropirate samples from mainbuff, 
                    //preparing to create spectrum
                    for (n = 0; n < len; n++)
                    {
                        DFT_inreal[n] = origdata[(int)(n + from)];
                    }

                    // non 2^n spectrum
                    dft(DFT_inreal, DFT_inimag, DFT_outreal, DFT_outimag);

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 0; n < 51 * 2; n += 2)
                    {
                        intpdata_big[n] = DFT_outreal[n / 2] / comp;
                        intpdata_big[n + 1] = DFT_outimag[n / 2] / comp;
                    }

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 53 * 2 - 1; n < DFT_outreal.Length * 2; n += 2)
                    {
                        intpdata_big[n + 921 * 2 + 1] = DFT_outreal[n / 2] / comp;
                        intpdata_big[n + 1 + 921 * 2 + 1] = DFT_outimag[n / 2] / comp;
                    }

                    // IFFT of widened spectrum
                    fft(intpdata_big, 1024, -1);

                    for (n = 0; n < 533; n++) //copy convert and saturate
                    {
                        if (intpdata_big[(n + 512 - 266) * 2] > 255.0)
                            intpdata[n] = 255;
                        else if (intpdata_big[(n + 512 - 266) * 2] < 0.0)
                            intpdata[n] = 0;
                        else
                            intpdata[n] = (byte)intpdata_big[(n + 512 - 266) * 2];

                    }

                    break;

                case 20:   //20x

                    // in final 533 points are needed. Nearest bigger 2^n is 1024
                    // so 1024/20 = 51.2
                    len = 51;
                    comp = 51.2;
                    from = 266.0 - len / 2.0;

                    DFT_inreal = new double[len];
                    DFT_inimag = new double[len];
                    DFT_outreal = new double[len];
                    DFT_outimag = new double[len];


                    // taking appropirate samples from mainbuff, 
                    //preparing to create spectrum
                    for (n = 0; n < len; n++)
                    {
                        DFT_inreal[n] = origdata[(int)(n + from)];
                    }

                    // non 2^n spectrum
                    dft(DFT_inreal, DFT_inimag, DFT_outreal, DFT_outimag);

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 0; n < 24 * 2 + 1; n += 2)
                    {
                        intpdata_big[n] = DFT_outreal[n / 2] / comp;
                        intpdata_big[n + 1] = DFT_outimag[n / 2] / comp;
                    }

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 27 * 2; n < DFT_outreal.Length * 2; n += 2)
                    {
                        intpdata_big[n + 973 * 2] = DFT_outreal[n / 2] / comp;
                        intpdata_big[n + 1 + 973 * 2] = DFT_outimag[n / 2] / comp;
                    }

                    // IFFT of widened spectrum
                    fft(intpdata_big, 1024, -1);

                    for (n = 0; n < 533; n++) //copy convert and saturate
                    {
                        if (intpdata_big[(n + 512 - 266) * 2] > 255.0)
                            intpdata[n] = 255;
                        else if (intpdata_big[(n + 512 - 266) * 2] < 0.0)
                            intpdata[n] = 0;
                        else
                            intpdata[n] = (byte)intpdata_big[(n + 512 - 266) * 2];

                    }

                    break;


                case 50:   //50x

                    // in final 533 points are needed. Nearest bigger 2^n is 1024
                    // so 1024/50 = 20.48
                    len = 20;
                    comp = 20.48;
                    from = 266.0 - len / 2.0;

                    DFT_inreal = new double[len];
                    DFT_inimag = new double[len];
                    DFT_outreal = new double[len];
                    DFT_outimag = new double[len];


                    // taking appropirate samples from mainbuff, 
                    //preparing to create spectrum
                    for (n = 0; n < len; n++)
                    {
                        DFT_inreal[n] = origdata[(int)(n + from)];
                    }

                    // non 2^n spectrum
                    dft(DFT_inreal, DFT_inimag, DFT_outreal, DFT_outimag);

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 0; n < 10 * 2; n += 2)
                    {
                        intpdata_big[n] = DFT_outreal[n / 2] / comp;
                        intpdata_big[n + 1] = DFT_outimag[n / 2] / comp;
                    }

                    //creating the spectrum widening
                    //lower half
                    // and compensating the amplitude
                    for (n = 12 * 2 - 1; n < DFT_outreal.Length * 2; n += 2)
                    {
                        intpdata_big[n + 1003 * 2 + 1] = DFT_outreal[n / 2] / comp;
                        intpdata_big[n + 1 + 1003 * 2 + 1] = DFT_outimag[n / 2] / comp;
                    }

                    // IFFT of widened spectrum
                    fft(intpdata_big, 1024, -1);

                    for (n = 0; n < 533; n++) //copy convert and saturate
                    {
                        if (intpdata_big[(n + 512 - 266) * 2] > 255.0)
                            intpdata[n] = 255;
                        else if (intpdata_big[(n + 512 - 266) * 2] < 0.0)
                            intpdata[n] = 0;
                        else
                            intpdata[n] = (byte)intpdata_big[(n + 512 - 266) * 2];

                    }

                    break;

            }

        }
    }
}