using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

namespace vcsCALFSR
{
    public partial class Form1 : Form
    {
        public Form1()
        {
            InitializeComponent();
        }
        //--------------------------------------------------------------------

        
        private void Form1_Load(object sender, EventArgs e)
        {
            toolStripLabel1.Text = "バンド幅";
            toolStripTextBox1.Text = "0.4";
        }
        //--------------------------------------------------------------------

        
        private void toolStripButton1_Click(object sender, EventArgs e)
        {
            main_part();
            MessageBox.Show("処理完了", "完了通知");
        }
        //--------------------------------------------------------------------
        
        
        private void main_part()
        {
            //フーリエスペクトル比算定
            //分子・分母にするデータのサンプリングピッチ・データ数は等しくとる
            int nnn,i;
            string fnameR = "";
            string fnameW = "";
            System.IO.StreamReader sr;
            System.IO.StreamWriter sw;
            string dat="";
            string[] sbuf;
            char delim = ',';
            int ncase,icase;
            double dt=0.0,band;
            int nn=0,nd;

            int ndmax = (int)Math.Pow(2,16);
            double[] fs01=new double[ndmax+1];
            double[] fs02=new double[ndmax+1];
            double[] fsp1=new double[ndmax+1];
            double[] fsp2=new double[ndmax+1];
            double[] fq=new double[ndmax+1];

            for(i=0;i<=ndmax;i++){
                fs01[i] = 0.0;
                fs02[i] = 0.0;
                fsp1[i] = 0.0;
                fsp2[i] = 0.0;
                fq[i] = 0.0;
            }

            //平滑化バンド幅設定
            band = double.Parse(toolStripTextBox1.Text);

            //入出力ファイル読み込み
            openFileDialog1.InitialDirectory = System.IO.Directory.GetCurrentDirectory();
            if(openFileDialog1.ShowDialog() == DialogResult.OK)fnameR = openFileDialog1.FileName;
            sr = new System.IO.StreamReader(fnameR, System.Text.Encoding.Default);
            dat = sr.ReadLine() ; sbuf = dat.Split(delim) ; ncase = int.Parse(sbuf[0]);
            string[] strcom1=new string[ncase];
            string[] strcom2=new string[ncase];
            string[] rn1=new string[ncase];
            string[] rn2=new string[ncase];
            string[] wn0=new string[ncase];
            for(icase=1;icase<=ncase;icase++){
                dat = sr.ReadLine() ; sbuf = dat.Split(delim);
                strcom1[icase - 1] = sbuf[0]; //分子スペクトル名
                strcom2[icase - 1] = sbuf[1]; //分母スペクトル名
                rn1[icase - 1] = sbuf[2]; //分子加速度時刻歴ファイル名
                rn2[icase - 1] = sbuf[3]; //分母加速度時刻歴ファイル名
                wn0[icase - 1] = sbuf[4]; //結果出力ファイル名
            }
            sr.Close();

            for(icase=1;icase<=ncase;icase++){
                for(nnn=1;nnn<=2;nnn++){
                    //原波形読み込みと平滑化フーリエスペクトル計算
                    switch(nnn){
                        case 1: //分子側加速度時刻歴入力・スペクトル計算
                            fnameR = System.IO.Path.GetDirectoryName(fnameR) + "\\" + rn1[icase - 1];
                            CALFSP(fnameR, out dt, out nn, band, fs01, fsp1);
                            break;
                        case 2: //分母側加速度時刻歴入力・スペクトル計算
                            fnameR = System.IO.Path.GetDirectoryName(fnameR) + "\\" + rn2[icase - 1];
                            CALFSP(fnameR, out dt, out nn, band, fs02, fsp2);
                            break;
                    }
                }
                //スペクトル比計算・振動数設定
                nd = (int)(nn / 2);
                double[] fspr=new double[nd+1];
                for(i=0;i<=nd;i++){
                    fspr[i] = fsp1[i] / fsp2[i];
                    fq[i] = (double)i / (double)nn / dt; //振動数
                }

                //結果書き出し
                fnameW = System.IO.Path.GetDirectoryName(fnameR) + "\\" + wn0[icase - 1];
                sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
                dat = "分子ファイル名" + "," + strcom1[icase - 1] + "," + rn1[icase - 1] ; sw.WriteLine(dat);
                dat = "分母ファイル名" + "," + strcom2[icase - 1] + "," + rn2[icase - 1] ; sw.WriteLine(dat);
                dat = "出力ファイル名" + "," + "sp比" + "," + wn0[icase - 1] ; sw.WriteLine(dat);
                dat = "バンド幅" + "," + band.ToString() ; sw.WriteLine(dat);
                dat = "データ数" + "," + (nd + 1).ToString() ; sw.WriteLine(dat);
                dat = "振動数,分子原SP,分母原SP,平滑化分子SP,平滑化分母SP,平滑化SP比" ; sw.WriteLine(dat);
                for(i=0;i<=nd;i++){
                    dat = fq[i].ToString("E");
                    dat = dat + "," + fs01[i].ToString("E");
                    dat = dat + "," + fs02[i].ToString("E");
                    dat = dat + "," + fsp1[i].ToString("E");
                    dat = dat + "," + fsp2[i].ToString("E");
                    dat = dat + "," + fspr[i].ToString("E");
                    sw.WriteLine(dat);
                }
                sw.Close();
            }

            //出力ファイル一覧出力
            fnameW = System.IO.Path.GetDirectoryName(fnameR) + "\\fileoutlist.txt";
            sw = new System.IO.StreamWriter(fnameW, false, System.Text.Encoding.Default);
            dat = ncase.ToString("0") ; sw.WriteLine(dat);
            for(icase=1;icase<=ncase;icase++){
                dat = wn0[icase - 1] ; sw.WriteLine(dat);
            }
            sw.Close();

        }

        //--------------------------------------------------------------------

        
        private void CALFSP(string fnameR, out double dt, out int nn, double band, double[] fs0, double[] fsp)
        {
            System.IO.StreamReader sr;
            string dat;
            string[] sbuf;
            char delim = ',';
            int i, ndata, nd;
            double ddymax, dymax, ymax, df;

            //原波形読み込み
            sr = new System.IO.StreamReader(fnameR, System.Text.Encoding.Default);
            dat = sr.ReadLine(); sbuf = dat.Split(delim);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); dt = double.Parse(sbuf[1]);
            dat = sr.ReadLine(); sbuf = dat.Split(delim); ndata = int.Parse(sbuf[1]);
            nn = 2;
            do
            {
                nn = nn * 2;
            } while (nn < ndata);
            nd = (int)(nn / 2);
            double[] ddy = new double[nn];
            double[] dy = new double[nn];
            double[] y = new double[nn];
            double[] xr = new double[nn];
            double[] xi = new double[nn];
            for (i = 1; i <= nn; i++)
            { //データ数を2の階乗にセット
                ddy[i - 1] = 0.0; dy[i - 1] = 0.0; y[i - 1] = 0.0;
            }
            for (i = 1; i <= ndata; i++)
            {
                dat = sr.ReadLine(); sbuf = dat.Split(delim); ddy[i - 1] = double.Parse(sbuf[0]);
            }
            sr.Close();

            //平滑化フーリエスペクトル作成
            IACC(ndata, dt, ddy, dy, y, out ddymax, out dymax, out ymax); //数値積分
            CRAC(ndata, dt, ddy, dy, y, ddymax, dymax, ymax); //基線補正
            for (i = 1; i <= nn; i++)
            {
                xr[i - 1] = ddy[i - 1];
                xi[i - 1] = 0.0;
            }
            FFT(nn, xr, xi); //フーリエ変換
            for (i = 1; i <= nn; i++)
            { //戻り値をデータ数で除す
                xr[i - 1] = xr[i - 1] / (double)nn;
                xi[i - 1] = xi[i - 1] / (double)nn;
            }
            for (i = 1; i <= nd + 1; i++)
            {
                fs0[i - 1] = dt * (double)nn * Math.Sqrt(xr[i - 1] * xr[i - 1] + xi[i - 1] * xi[i - 1]); //分子フーリエスペクトル
                fsp[i - 1] = fs0[i - 1];
            }
            df = 1.0 / dt / (double)nn;
            SWIN(nn, fsp, df, band); //スペクトル平滑化

        }
        //------------------------------------------------------------------------------

        
        private void FFT(int nn, double[] xr, double[] xi)
        {
            //**************************************
            //高速フーリエ変換・逆変換
            //**************************************
            //nn  ：データ数（2の累乗）
            //xr[]：入出力データ実数部
            //xi[]：入出力データ虚数部

            int g, h, i, j, k, l, m, n, p, q;
            double a, b, xd;

            n = nn;
            double[] s = new double[(int)(n / 2) + 1];
            double[] c = new double[(int)(n / 2) + 1];

            i = 0; j = 0; k = 0; l = 0; p = 0; h = 0; g = 0; q = 0;
            m = (int)(Math.Log((double)n) / Math.Log(2.0)+1.0);
            a = 0.0;
            b = Math.PI * 2.0 / (double)n;
            for (i = 0; i <= (int)(n / 2); i++)
            {
                s[i] = Math.Sin(a);
                c[i] = Math.Cos(a);
                a = a + b;
            }
            l = n; h = 1;
            for (g = 1; g <= m; g++)
            {
                l = (int)(l / 2); k = 0;
                for (q = 1; q <= h; q++)
                {
                    p = 0;
                    for (i = k; i <= l + k - 1; i++)
                    {
                        j = i + l;
                        a = xr[i] - xr[j];
                        b = xi[i] - xi[j];
                        xr[i] = xr[i] + xr[j];
                        xi[i] = xi[i] + xi[j];
                        if (p == 0)
                        {
                            xr[j] = a;
                            xi[j] = b;
                        }
                        else
                        {
                            xr[j] = a * c[p] + b * s[p];
                            xi[j] = b * c[p] - a * s[p];
                        }
                        p = p + h;
                    }
                    k = k + l + l;
                }
                h = h + h;
            }
            j = (int)(n / 2);
            for (i = 1; i <= n - 1; i++)
            {
                k = n;
                if (j < i)
                {
                    xd = xr[i]; xr[i] = xr[j]; xr[j] = xd;
                    xd = xi[i]; xi[i] = xi[j]; xi[j] = xd;
                }
                k = (int)(k / 2);
                while (j >= k)
                {
                    j = j - k; k = (int)(k / 2);
                    if (k == 0) break;
                }
                j = j + k;
            }
        }
        //---------------------------------------------------------------------------------

        private void IACC(int nn, double dt, double[] ddy, double[] dy, double[] y,
                          out double ddymax, out double dymax, out double ymax)
        {
            //**************************************
            //加速度時刻歴数値積分
            //**************************************
            //nn:加速度時刻歴データ総数
            //dt    :時間間隔（入力値）
            //ddy[] :加速度時刻歴（入力値）
            //dy[]  :速度時刻歴（計算出力）
            //y[]   :変位時刻歴（計算出力値）
            //ddymax:加速度最大値（計算出力値）
            //dymax :速度最大値（計算出力値）
            //ymax  :変位最大値（計算出力値）

            int i;

            ddymax = 0.0;
            dymax = 0.0;
            ymax = 0.0;
            dy[0] = 0.0;
            y[0] = 0.0;
            for (i = 1; i <= nn - 1; i++)
            {
                dy[i] = dy[i - 1] + (ddy[i - 1] + ddy[i]) * dt / 2.0;
                y[i] = y[i - 1] + dy[i - 1] * dt + (ddy[i - 1] / 3.0 + ddy[i] / 6.0) * dt * dt;
                if (ddymax < Math.Abs(ddy[i])) ddymax = Math.Abs(ddy[i]);
                if (dymax < Math.Abs(dy[i])) dymax = Math.Abs(dy[i]);
                if (ymax < Math.Abs(y[i])) ymax = Math.Abs(y[i]);
            }
        }
        //---------------------------------------------------------------------------
        private void CRAC(int nn, double dt, double[] ddy, double[] dy, double[] y,
                          double ddymax, double dymax, double ymax)
        {
            //**************************************
            //加速度時刻歴基線補正
            //**************************************
            //nn:加速度時刻歴データ総数
            //dt    :時間間隔（入力値）
            //ddy[] :加速度時刻歴（入力値・書換出力値）
            //dy[]  :速度時刻歴（入力値）
            //y[]   :変位時刻歴（入力値）
            //ddymax:加速度最大値（入力値）
            //dymax :速度最大値（入力値）
            //ymax  :変位最大値（入力値）

            int i;
            double tt, t, sum, a1, a0, acmax, coef;

            tt = (double)(nn - 1) * dt;
            t = 0.0;
            for (i = 0; i <= nn - 1; i++)
            {
                y[i] = y[i] * (3.0 * tt - 2.0 * t) * t * t;
                t = t + dt;
            }
            sum = (y[0] + y[nn - 1]) / 2.0;
            for (i = 1; i <= nn - 2; i++)
            {
                sum = sum + y[i];
            }
            sum = sum * dt;
            a1 = 28.0 / 13.0 / tt / tt * (2.0 * dy[nn - 1] - 15.0 / Math.Pow(tt, 5.0) * sum);
            a0 = dy[nn - 1] / tt - a1 / 2.0 * tt;
            t = 0.0;
            acmax = 0.0;
            for (i = 0; i <= nn - 1; i++)
            {
                ddy[i] = ddy[i] - a0 - a1 * t;
                if (acmax < Math.Abs(ddy[i])) acmax = Math.Abs(ddy[i]);
                t = t + dt;
            }
            coef = ddymax / acmax;
            for (i = 0; i <= nn - 1; i++)
            {
                ddy[i] = ddy[i] * coef;
            }
        }
        //---------------------------------------------------------------------------

        private void SWIN(int nn, double[] fs, double df, double band)
        {
            //**************************************
            //スペクトル平滑化
            //**************************************
            //nn   ：時刻歴全数
            //fs() ：フーリエ・スペクトル
            //df   ：スペクトルの振動数間隔（Hz）
            //band ：バンド幅（Hz）
            //nfold：スペクトル値総数

            double pi = Math.PI;
            double tt, udf, dif, s;
            int i, j;
            int nfold, lmax, ll, ln, lt, le;
            double[] w = new double[nn];
            double[] g = new double[nn];
            double[] g1 = new double[nn];
            double[] g2 = new double[nn];

            nfold = (int)(nn / 2) + 1;

            if (band > 0.0)
            {
                tt = 1.0 / df;
                udf = 280.0 / 151.0 / band * df;
                lmax = (int)(2.0 / udf) + 1;

                if (udf <= 0.5)
                {
                    w[0] = 0.75 * udf;
                    for (i = 2; i <= lmax; i++)
                    {
                        dif = 0.5 * pi * (double)(i - 1) * udf;
                        w[i - 1] = w[0] * Math.Pow(Math.Sin(dif) / dif, 4.0);
                    }
                    g[0] = fs[0] * fs[0] / tt;
                    for (i = 2; i <= nfold - 1; i++)
                    {
                        g[i - 1] = 2.0 * fs[i - 1] * fs[i - 1] / tt;
                    }
                    g[nfold - 1] = fs[nfold - 1] * fs[nfold - 1] / tt;
                    ll = lmax * 2 - 1;
                    ln = ll - 1 + nfold;
                    lt = (ll - 1) * 2 + nfold;
                    le = lt - lmax + 1;
                    for (i = 1; i <= lt; i++)
                    {
                        g1[i - 1] = 0.0;
                    }
                    for (i = 1; i <= nfold; i++)
                    {
                        g1[ll - 2 + i] = g[i - 1];
                    }
                    for (i = lmax; i <= le; i++)
                    {
                        s = w[0] * g1[i - 1];
                        for (j = 2; j <= lmax; j++)
                        {
                            s = s + w[j - 1] * (g1[i - j] + g1[i + j - 2]);
                        }
                        g2[i - 1] = s;
                    }
                    for (j = 2; j <= lmax; j++)
                    {
                        g2[ll + j - 2] = g2[ll + j - 2] + g2[ll - j];
                        g2[ln - j] = g2[ln - j] + g2[ln + j - 2];
                    }
                    for (i = 1; i <= nfold; i++)
                    {
                        g[i - 1] = g2[ll - 2 + i];
                    }
                    fs[0] = Math.Sqrt(g[0] * tt);
                    for (i = 2; i <= nfold - 1; i++)
                    {
                        fs[i - 1] = Math.Sqrt(g[i - 1] * tt / 2.0);
                    }
                    fs[nfold - 1] = Math.Sqrt(g[nfold - 1] * tt);
                }
            }
        }
        //--------------------------------------------------------------------

    }
}