//====================================================================
#include <vcl.h>
#pragma hdrstop
//debug 2005/03/12 Hal.T

#include "sindo.h"
//====================================================================
Sindo_Base::Sindo_Base()
{
  fft = 0;
  Kfil  = 0;
  buf = 0;
  mbuf  = 0;
  vbuf  = 0;
  abuf  = 0;    //debug 2005/03/12 Hal.T
  Init();
}
//----------------------------------------------
Sindo_Base::~Sindo_Base()
{
  if(fft ) delete  fft;
  if(Kfil) delete  []  Kfil;
  if(buf ) delete  []  buf;
  if(mbuf) delete  []  mbuf;
  if(vbuf) delete  []  vbuf;
  if(abuf) delete  []  abuf;

}
//----------------------------------------------
//初期化
void  Sindo_Base::Init(void)
{
  LN  = 11;   //FFT bit  これを変化させると計算サイズが内部変化する

  AN  = pow(2,LN); //2^LN  //1024
  dT  = 0.01; //毎秒１００サンプル
  interval  = 50; //計測震度の計算間隔、サンプル数単位で表現する
  it  = 0;        //計測震度の計算勘定

  //--------------------
  //バッファの初期化
  if(buf) delete  []  buf;
  buf = new V3D  [AN];

  if(vbuf) delete  []  vbuf;
  vbuf = new V3D  [AN];

  if(abuf) delete  []  abuf;
  abuf = new double [AN];

  if(mbuf) delete  []  mbuf;
  mbuf = new V3D  [interval];
  mbufN = 0;


  //--------------------
  //FFTの初期化
  if(fft) delete  fft;
  fft = new FFT(LN);
  KeisokuSindo  = 0.0;
  KCal = 0;

  //--------------------
  //フィルタの初期化
  if(Kfil) delete  []  Kfil;
  Kfil = new double  [AN];
  double  FreqStep =  1.0/(dT*AN);
  int i;
  double  f, d;

  Kfil[0] = 0.0;

  for(i = 1;  i<AN/2; i++)  {
    f = FreqStep * i;
    d = sqrt(1/f)*pow(1+0.694*pow(f/10,2)+0.241*pow(f/10,4)+0.0557*pow(f/10,6)
        +0.009664*pow(f/10,8)+0.00134*pow(f/10,10)+0.000155*pow(f/10,12),-1/2)
        *pow(1-exp(-pow(f/(0.5),3)),0.5);
    Kfil[i] = d;
  }

  Kfil[i] = 0.0;  //SPS/2

  for(i=1 ;  i<AN/2; i++)  {
    Kfil[AN/2+i] = Kfil[AN/2-i];
  }
}
//----------------------------------------------
//データ投入　単位はgal
//震度が更新されたら 1、無更新なら０で帰る
int Sindo_Base::set_gal(V3D  gal)
{
  return  set_gal(&gal, 1);
}
//----------------------------------------------
/*データ投入　単位はgal
  震度が更新されたら 1、無更新なら０で帰る
  未計算バッファに突っ込んで、いっぱいになるたびに計算する
  計測震度値はこの間で最大のものを維持する
*/
int Sindo_Base::set_gal(V3D *gal,  int n)
{
  double  KST;  //計測震度の一時記憶
  KCal = 0;
  KST = KeisokuSindo  = 0.0;

  while(n>0)  {
    //一時バッファへのコピー
    //コピー可能量 interval - mbufN
    int mcv = interval - mbufN;
    int cpN;  //コピーする個数

    if(n <= mcv)  {
      cpN = n;
    } else  {
      cpN = mcv;
    }
    memcpy(&mbuf[mbufN], gal, cpN*sizeof(V3D));
    mbufN +=  cpN;
    gal   +=  cpN;  //次へ進めておく
    n     -=  cpN;

    if(mbufN  ==  interval) { //計測震度の計算をする
      //過去データの移動
      memmove(&buf[0],  &buf[interval], (AN - interval)*sizeof(V3D));
      //一時バッファからのコピー
      memcpy(&buf[AN - interval], mbuf, interval*sizeof(V3D));
      mbufN = 0;  //一時バッファの数をクリアしておく

      //計測震度の計算
      KST = Calc();   //計測震度の計算
      if(KST>KeisokuSindo)  KeisokuSindo = KST;
      KCal = 1;
    }
  }
  return  KCal;
}
//----------------------------------------------
//計測震度を得る、計算そのものは終わっている
double  Sindo_Base::get_sindo(void)
{
  return  KeisokuSindo;
}
//----------------------------------------------
//計算の実行
double  Sindo_Base::Calc(void)
{
  //buffからFFTに送ってフィルタ計算をする
  //X軸
  for(int n=0;  n<AN; n++)  {
    fft->src[n].x = buf[n].x;    fft->src[n].y = 0.0;
  }
  fft->fft_cnv(1);   //順FFT
  //マスクをかける
  for(int n=0;  n<AN; n++)  {
    fft->src[n] *= Kfil[n];   //複素×ダブル
  }
  fft->fft_cnv(-1);   //逆FFT

  for(int n=0;  n<AN; n++)  {
    vbuf[n].x = fft->src[n].abs();    //複素になっている、絶対値として扱おう
  }

  //Y軸
  for(int n=0;  n<AN; n++)  {
    fft->src[n].x = buf[n].y;    fft->src[n].y = 0.0;
  }
  fft->fft_cnv(1);   //順FFT
  //マスクをかける
  for(int n=0;  n<AN; n++)  {
    fft->src[n] *= Kfil[n];   //複素×ダブル
  }
  fft->fft_cnv(-1);   //逆FFT

  for(int n=0;  n<AN; n++)  {
    vbuf[n].y = fft->src[n].abs();    //複素になっている、絶対値として扱おう
  }

  //Z軸
  for(int n=0;  n<AN; n++)  {
    fft->src[n].x = buf[n].z;    fft->src[n].y = 0.0;
  }
  fft->fft_cnv(1);   //順FFT
  //マスクをかける
  for(int n=0;  n<AN; n++)  {
    fft->src[n] *= Kfil[n];   //複素×ダブル
  }
  fft->fft_cnv(-1);   //逆FFT

  for(int n=0;  n<AN; n++)  {
    vbuf[n].z = fft->src[n].abs();    //複素だ、絶対値として扱おう
  }

  //ベクトル合成と規格化
  double  kikaku = 1.0/AN;
  for(int n=0;  n<AN; n++)  abuf[n] = vbuf[n].abs() * kikaku;

  /*継続時間を考慮した振幅の決定
    最大値を調べて、最大値をてがかりにヒストグラムを作成

  */
  double  tansaku, smax = 0.0;
  for(int n=0;  n<AN; n++)  if(smax < abuf[n])  smax = abuf[n];

  if(smax==0.0) return  KeisokuSindo = 0.0;

  int tn300 = taw/dT;   //継続時間分のカウント数
  const int divh = 4000;
  int histgram[divh+1];   //分解能４０００でヒストグラムを作る

  for(int n=0;  n<divh+1; n++)  histgram[n]=0;

  for(int n=0;  n<AN; n++)  {
    histgram[(int)(((double)divh * abuf[n])/smax)]++;
  }
  int n;
  for(n=divh-1;  n>=0; n--)  {
    histgram[n]+=histgram[n+1];
    if(histgram[n] >= tn300)  break;
  }

  double  kg = n * smax / divh;
  KeisokuSindo = 2.0 * log10(kg) + 0.94;
  if(KeisokuSindo < 0.0) KeisokuSindo = 0.0;

  KeisokuSindo = round(KeisokuSindo, 1);  //vecmath.h

  return  KeisokuSindo;
}
//====================================================================
static  char  *IniFileName = "NoiseData.csv";

Sindo::Sindo() : Sindo_Base()
{
  d_num = noise_cal = 0;
  NoiseData = new V3D [AN];  //ポインタの初期化
  for(int n=0;  n<AN; NoiseData[n++]=V3D(0.0, 0.0, 0.0));

  //Make File name
  //カレントディレクトリのフルパスを作成する
  GetCurrentDirectory(MAXPATH,  fname);
  if(fname[strlen(fname)-1] != '\\')  strcat(fname,  "\\");
  strcat(fname, IniFileName);

  FILE  *fp;
  fp = fopen("NoiseData.csv", "r");
  if(fp)  {
    fclose(fp);
    ReadCSV("NoiseData.csv", NoiseData, AN);
  }
}
//----------------------------------------------
Sindo::~Sindo()
{
  WriteCSV(fname, NoiseData, AN);

  delete  []  NoiseData;
}
//----------------------------------------------
int Sindo::set_gal(V3D *gal,  int n)
{
//  noise_cal = 0;
  d_num +=  n;
  if(d_num > AN)  {
    d_num = 0;  noise_cal = 1;
  }
  return  Sindo_Base::set_gal(gal, n);
}
//----------------------------------------------
//計算の実行 新しいバージョンに置き換える　（Sindo_Base::Calc　はvirtual指示）
double  Sindo::Calc(void)
{
  V2D v2zero;
  //buffからFFTに送ってフィルタ計算をする
  //X軸  =============================
  for(int n=0;  n<AN; n++)  {
    fft->src[n].x = buf[n].x;    fft->src[n].y = 0.0;
  }
  fft->fft_cnv(1);   //順FFT

  if(noise_cal) {
    for(int n=1;  n<AN; n++)  {
      NoiseData[n].x = (9.0 * NoiseData[n].x + fft->src[n].abs())*0.1;
    }
  }
  for(int n=1;  n<AN; n++)  {
    double  vk = fft->src[n].abs();
    if(NoiseData[n].x < vk) {
      fft->src[n] *=  (vk*vk-NoiseData[n].x*NoiseData[n].x)/(vk*vk);
    } else  fft->src[n] = v2zero;
  }

  //マスクをかける
  for(int n=0;  n<AN; n++)  {
    fft->src[n] *= Kfil[n];   //複素×ダブル
  }
  fft->fft_cnv(-1);   //逆FFT

  for(int n=0;  n<AN; n++)  {
    vbuf[n].x = fft->src[n].abs();    //複素になっている、絶対値として扱おう
  }

  //Y軸  =============================
  for(int n=0;  n<AN; n++)  {
    fft->src[n].x = buf[n].y;    fft->src[n].y = 0.0;
  }
  fft->fft_cnv(1);   //順FFT

  if(noise_cal) {
    for(int n=1;  n<AN; n++)  {
      NoiseData[n].y = (9.0 * NoiseData[n].y + fft->src[n].abs())*0.1;
    }
  }
  for(int n=1;  n<AN; n++)  {
    double  vk = fft->src[n].abs();
    if(NoiseData[n].y < vk) {
      fft->src[n] *=  (vk*vk-NoiseData[n].y*NoiseData[n].y)/(vk*vk);
    } else  fft->src[n] = v2zero;
  }

  //マスクをかける
  for(int n=0;  n<AN; n++)  {
    fft->src[n] *= Kfil[n];   //複素×ダブル
  }
  fft->fft_cnv(-1);   //逆FFT

  for(int n=0;  n<AN; n++)  {
    vbuf[n].y = fft->src[n].abs();    //複素になっている、絶対値として扱おう
  }

  //Z軸  =============================
  for(int n=0;  n<AN; n++)  {
    fft->src[n].x = buf[n].z;    fft->src[n].y = 0.0;
  }
  fft->fft_cnv(1);   //順FFT

  if(noise_cal) {
    for(int n=1;  n<AN; n++)  {
      NoiseData[n].z = (9.0 * NoiseData[n].z + fft->src[n].abs())*0.1;
    }
  }
  for(int n=1;  n<AN; n++)  {
    double  vk = fft->src[n].abs();
    if(NoiseData[n].z < vk) {
      fft->src[n] *=  (vk*vk-NoiseData[n].z*NoiseData[n].z)/(vk*vk);
    } else  fft->src[n] = v2zero;
  }


  //マスクをかける
  for(int n=0;  n<AN; n++)  {
    fft->src[n] *= Kfil[n];   //複素×ダブル
  }
  fft->fft_cnv(-1);   //逆FFT

  for(int n=0;  n<AN; n++)  {
    vbuf[n].z = fft->src[n].abs();    //複素だ、絶対値として扱おう
  }

  //ベクトル合成と規格化
  double  kikaku = 1.0/AN;
  for(int n=0;  n<AN; n++)  abuf[n] = vbuf[n].abs() * kikaku;

  /*継続時間を考慮した振幅の決定
    最大値を調べて、最大値をてがかりにヒストグラムを作成

  */
  double  tansaku, smax = 0.0;
  for(int n=0;  n<AN; n++)  if(smax < abuf[n])  smax = abuf[n];

  if(smax==0.0) return  KeisokuSindo = 0.0;

  int tn300 = taw/dT;   //継続時間分のカウント数
  const int divh = 4000;
  int histgram[divh+1];   //分解能４０００でヒストグラムを作る

  for(int n=0;  n<divh+1; n++)  histgram[n]=0;

  for(int n=0;  n<AN; n++)  {
    histgram[(int)(((double)divh * abuf[n])/smax)]++;
  }
  int n;
  for(n=divh-1;  n>=0; n--)  {
    histgram[n]+=histgram[n+1];
    if(histgram[n] >= tn300)  break;
  }

  double  kg = n * smax / divh;
  KeisokuSindo = 2.0 * log10(kg) + 0.94;
  if(KeisokuSindo < 0.0) KeisokuSindo = 0.0;

  KeisokuSindo = round(KeisokuSindo, 2);  //vecmath.h　＠

  noise_cal = 0;

  return  KeisokuSindo;
}
//----------------------------------------------
void  Sindo::reset_noise(void)
{
  V3D zero; // = all zero
  for(int n=1;  n<AN; n++) NoiseData[n] = zero;
}


//====================================================================


#pragma package(smart_init)