//---------------------------------------------------------------------------
#include <vcl.h>
//#include <stdio.h>
//#include <math.h>
#include <complex.h>
#include <except.h>
#pragma hdrstop

#include "main.h"
#include "about.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma link "cspin"
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
    : TForm(Owner)
{
}
//---------------------------------------------------------------------------
// 入力波形の選択

void __fastcall TForm1::Btn_input1Click(TObject *Sender)
{
    if ( OpenDialog1->Execute() ){
        Edit_input1->Text = OpenDialog1->FileName;
    }
}
//---------------------------------------------------------------------------
void __fastcall TForm1::FormCreate(TObject *Sender)
{
        ClientHeight = 215;
        ClientWidth  = 470;
}
//---------------------------------------------------------------------------
// 高速フーリエ変換
//      参考文献：大崎順彦「新・地震動のスペクトル解析入門」鹿島出版会

void TForm1::fast ( complex<double> x[], int nn, int ind )
{
complex<double> temp, theta;
int i, j, k, m, kmax, istep;
const double PI = 3.14159265358979 ;

    j = 0 ;
    for ( i=0; i<nn; i++ ) {
        if ( i<j ) {
            temp = x[j] ;
            x[j] = x[i] ;
            x[i] = temp ;
        }
        m = nn / 2 ;
        while ( j+1>m ) {
            j = j - m ;
            m = m / 2 ;
            if ( m<2 ) break;
        }
        j = j + m ;
    }
    kmax = 1 ;
    while ( kmax<nn ) {
        istep = kmax * 2 ;
        for ( k=0; k<kmax; k++ ) {
            theta = complex<double>( 0., PI * ind * k / (double)kmax ) ;
            for ( i=k; i<nn; i=i+istep ) {
                j = i + kmax ;
                temp = x[j] * exp( theta ) ;
                x[j] = x[i] - temp ;
                x[i] = x[i] + temp ;
            }
        }
        kmax = istep ;
    }
}
//---------------------------------------------------------------------------
// パラメータの入力忘れチェック

int TForm1::check_1()
{
int ok ;

    ok = 0 ;
    if ( Edit_dt->Text == "" ) { ok = 1 ;}
    if ( Edit_number->Text == "" ) { ok = 1 ;}
    if ( Edit_input1->Text == "" ) { ok = 1 ;}
    if ( Edit_input2->Text == "" ) { ok = 1 ;}
    if ( Edit_output->Text == "" ) { ok = 1 ;}
    if ( ok == 1 ) {
        error_1() ;
    }
    return ( ok ) ;
}
//---------------------------------------------------------------------------
// パラメータの入力忘れに対するメッセージボックス表示

void TForm1::error_1()
{
    Application->MessageBox("パラメーターが入力されていない箇所があります。",
        "エラー",MB_ICONEXCLAMATION | MB_OK ) ;
}
//---------------------------------------------------------------------------
void __fastcall TForm1::Btn_execClick(TObject *Sender)
{
float z;
double dt, a, *coh2, *freq, *xin, *pc_r, *pc_i, *time,
       *px_raw, *py_raw, *pc_raw_r, *pc_raw_i, *px, *py, *phase,
       tn, tn2, df ;
complex<double> *pc, *pc_raw, *pxx, *pyy ;
const complex<double> zero(0.,0.);
const double PI = 3.14159265358979 ;
int n, total, i, nn, nfold ;

    Form2->Hide() ;

// パラメーターが全部入力されているか
    if ( check_1() == 1 ) { return ; }

    dt  = StrToFloat( Edit_dt->Text );      // 時間刻み
    total = StrToInt( Edit_number->Text );  // 地震波個数

    a = log( (double)(total-1) ) / log( 2. ) + 1.0 ;
    nn = pow ( 2., (int) a ) ;
    nfold = nn / 2 ;
    tn = dt / (double)nn ;
    tn2 = 2. * tn ;
    df = 1. / (double)nn / dt ;

    try {                      // 例外のテスト
      xin = new double[nn];
      time = new double[nn];
      pxx = new complex<double>[nn];
      pyy = new complex<double>[nn];
      pc_raw = new complex<double>[nn];
      px_raw = new double[nn];
      py_raw = new double[nn];
      pc_raw_r = new double[nn];
      pc_raw_i = new double[nn];

    }
    catch (bad_alloc) { // bad_alloc が送出された場合にのみこのブロックに入る
        Application->MessageBox
        ("メモリが足りません。データ個数を減らして下さい。",
         "エラー",MB_ICONEXCLAMATION | MB_OK ) ;
        exit (-1);
    }

    for ( i=0; i<total; i++ )
        time[i] = dt * (double)i ;

// 波形１

    if (( fpin=fopen(Edit_input1->Text.c_str(),"rt")) == NULL ){
        Application->MessageBox("波形ファイルを設定し直して下さい。",
        "入力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
        return ;
    }

    i = 0;
    while( fscanf(fpin,"%f",&z ) != EOF ){
        xin[i] = (double) z;
        i++ ;
        if ( i >= total ) break;
    }
    fclose(fpin);
    n = i;
    senkai ( time, xin, n ) ;
    for ( i=n; i<nn; i++ ) xin[i] = 0. ;

    for ( i=0; i<nn; i++ )
        pxx[i] = complex<double>( xin[i], 0. ) ;
    fast  ( pxx, nn, -1 ) ;
    for ( i=0; i<nn; i++ ) {
        pxx[i] /= (double) nn ;
        px_raw[i] = pow( abs( pxx[i] ),2 ) * tn2 ;
    }

// 波形２

    if (( fpin=fopen(Edit_input2->Text.c_str(),"rt")) == NULL ){
        Application->MessageBox("波形ファイルを設定し直して下さい。",
        "入力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
        return ;
    }

    i = 0;
    while( fscanf(fpin,"%f",&z ) != EOF ){
        xin[i] = (double) z;
        i++ ;
        if ( i >= total ) break;
    }
    fclose(fpin);
    n = i;
    senkai ( time, xin, n ) ;
    for ( i=n; i<nn; i++ ) xin[i] = 0. ;

    for ( i=0; i<nn; i++ )
        pyy[i] = complex<double>( xin[i], 0. ) ;
    fast  ( pyy, nn, -1 ) ;
    for ( i=0; i<nfold; i++ ) {
        pyy[i] /= (double) nn ;
        py_raw[i] = pow( abs( pyy[i] ),2 ) * tn2 ;
    }

//  コヒーレンスとフェイズ

    for ( i=0; i<nfold; i++ ) {
        pc_raw[i] = pxx[i] * conj( pyy[i] ) * tn ;
        pc_raw_r[i] = real( pc_raw[i] ) ;
        pc_raw_i[i] = imag( pc_raw[i] ) ;
    }

    delete[] pxx ;
    delete[] pyy ;
    delete[] xin ;
    delete[] time ;

    px   = new double[nn];
    py   = new double[nn];
    pc_r = new double[nn];
    pc_i = new double[nn];

    filter( px_raw, px, nfold ) ;
    filter( py_raw, py, nfold ) ;
    filter( pc_raw_r, pc_r, nfold ) ;
    filter( pc_raw_i, pc_i, nfold ) ;

    delete[] px_raw ;
    delete[] py_raw ;
    delete[] pc_raw_r ;
    delete[] pc_raw_i ;

    freq = new double[nn];
    pc   = new complex<double>[nn];
    coh2 = new double[nn];
    phase= new double[nn];

    for ( i=0; i<nfold; i++ ) {
        freq[i] = df * (double)i ;
        pc[i] = complex<double> ( pc_r[i], pc_i[i] ) ;
        coh2[i] = pow( abs( pc[i] ), 2 ) / px[i] / py[i] * 4. ;
        phase[i] = atan( pc_i[i] / pc_r[i] ) * 180. / PI ;
        if ( pc_r[i]<0. ) {
          if ( pc_i[i] > 0. ) phase[i]=phase[i]+180.;
          else                phase[i]=phase[i]-180.;
        }
    }

// ファイルに出力

    if (( fpout=fopen(Edit_output->Text.c_str(),"wt")) == NULL ){
        Application->MessageBox("結果の出力ファイルを設定し直して下さい。",
        "出力ファイルエラー",MB_ICONEXCLAMATION | MB_OK ) ;
        return ;
    }

    fprintf(fpout,"%s\t%s\t%s\t%s\t%s\n",
        "周波数(Hz)","コヒーレンス","フェイズ","クロス（実数部）","クロス（虚数部）" );
    for ( i=0; i<nfold; i++ )
        fprintf(fpout,"%lf\t%lg\t%lg\t%lg\t%lg\n",
            freq[i], coh2[i], phase[i], pc_r[i], pc_i[i] );

    fclose(fpout);

    delete[] freq ;
    delete[] coh2 ;
    delete[] pc_r ;
    delete[] pc_i ;
    delete[] px ;
    delete[] py ;
    delete[] phase;

    Application->MessageBox("計算が終了しました。",
            "波形修正",MB_ICONEXCLAMATION | MB_OK ) ;

}
//---------------------------------------------------------------------------
// 三角フィルター
//                日野幹雄：スペクトル解析，朝倉書店，1977年．p.233

void TForm1::filter( double ssp[], double sp[], int m )
{
int k, i, k1, k2, nf = 10 ;
double sum ;

       for ( k=0; k<m; k++ ) {
           sum = ssp[k] * nf ;
           for ( i=1; i<nf; i++ ) {
               k1 = k - i ;
               k2 = k + i ;
               if ( k1 < 0 ) k1 = -k1 ;
               if ( k2 >= m ) k2 = 2 * m - k2 - 2 ;
               sum = sum + ( ssp[k1] + ssp[k2] ) * ( nf - i ) ;
           }
           sp[k] = sum / nf / nf ;
       }
}
//---------------------------------------------------------------------------
// 説明文の表示

void __fastcall TForm1::Btn_helpClick(TObject *Sender)
{
    Form2->Show() ;
}
//---------------------------------------------------------------------------
// 基線補正 線形回帰

void TForm1::senkai ( double x[], double y[], int n )
{
double ax=0., ay=0., b=0., c=0., alph, beta;
int i;

    for ( i=0; i<n; i++ ) {
        ax += x[i] ;
        ay += y[i] ;
    }
    ax /= n ;
    ay /= n ;

    for ( i=0; i<n; i++ ) {
        b += ( x[i] - ax ) * ( y[i] - ay ) ;
        c += ( x[i] - ax ) * ( x[i] - ax ) ;
    }
    beta = b / c ;
    alph = ay - beta * ax ;
    for ( i=0; i<n; i++ )
        y[i] = y[i] - alph - beta * x[i] ;
}
//---------------------------------------------------------------------------
// スペクトル保存

void __fastcall TForm1::Btn_specClick(TObject *Sender)
{
       if ( SaveDialog1->Execute() ){
        Edit_output->Text = SaveDialog1->FileName;
    }

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

void __fastcall TForm1::Btn_input2Click(TObject *Sender)
{
    if ( OpenDialog1->Execute() ){
        Edit_input2->Text = OpenDialog1->FileName;
    }

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