// 高速フーリエ変換
//      参考文献：大崎順彦「新・地震動のスペクトル解析入門」鹿島出版会

void TForm1::fast ( struct complex<double> x[], int nn, int ind )
{
struct 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 / 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 ;
	}
}
//---------------------------------------------------------------------------
// スペクトルの平滑化　Smoothing of Spectra by Parzen's Spectral Window
//      参考文献：大崎順彦「新・地震動のスペクトル解析入門」鹿島出版会

int TForm1::smooth ( int n, double f[], double g[],
    double df, double band, double g1[], double g2[] )
{
double dif, s, w[101] ;
int k, i, imax;
const double PI = 3.14159265358979 ;

// 初期化
	double t = 1.0 / df ;
	if ( band >= 0.0001 ) {
		double udf = 280. / 151. / band * df ;
		if ( udf > 0.5 ) {
			Application->MessageBox("バンド幅が狭すぎます。",
			"エラー",MB_ICONEXCLAMATION | MB_OK ) ;
			return ( -1 ) ;
		}
		imax = (int)( 2.0 / udf ) + 1 ;
		if ( imax > 101 ) {
			Application->MessageBox("バンド幅が広すぎます。",
			"エラー",MB_ICONEXCLAMATION | MB_OK ) ;
			return ( -1 ) ;
		}

// スペクトルウィンドウ
		w[0] = 3. / 4. * udf ;
		for ( i=1; i<imax; i++ ) {
			dif = PI / 2. * i * udf ;
			w[i] = 3. / 4. * udf * pow( sin(dif) / dif , 4 ) ;
		}
	}

// フーリエスペクトルからパワースペクトルを計算
	g[0] = pow( f[0], 2 ) / t ;
	for ( k=1; k<n; k++ )
		g[k] = 2.0 * pow( f[k], 2 ) / t ;
	g[n] = pow( f[n], 2 ) / t ;
	if ( band < 0.0001 ) return( 0 ) ;

// パワースペクトルのスムージング
	int ii = imax * 2 - 1 ;
	int in = ii + n ;
	int it = ( ii - 1 ) * 2 + n + 1 ;
	int ie = it - imax + 1 ;
	for ( k=0; k<it; k++ ) g1[k] = 0. ;
	for ( k=0; k<=n; k++ ) g1[ii-1+k] = g[k] ;
	for ( k=imax-1; k<ie; k++ ) {
		s = w[0] * g1[k] ;
		for ( i=1; i<imax; i++ )
			s = s + w[i] * ( g1[k-i] + g1[k+i] ) ;
		g2[k] = s ;
	}
	for ( i=1; i<imax; i++ ) {
		g2[ii+i-1] = g2[ii+i-1] + g2[ii-i-1] ;
		g2[in-i-1] = g2[in-i-1] + g2[in+i-1] ;
	}
	for ( k=0; k<=n; k++ ) g[k] = g2[ii-1+k] ;

// スムージングされたフーリエスペクトル
	f[0] = sqrt( g[0] * t ) ;
	for ( k=1; k<n; k++ )
	f[k] = sqrt( g[k] * t / 2.0 ) ;
	f[n] = sqrt( g[n] * t ) ;

	return ( 0 ) ;

}
//---------------------------------------------------------------------------
void main( )
{
float z;
double dt, a, band, df,
	*time, *freq, *accr, *temp1, *temp2,
	*afft, *pow_spec ;
struct complex<double> *comp ;
const struct complex<double> zero(0.,0.);
int n, total, i, nn, nfold ;

// ファイルのオープン

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

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

	dt = StrToFloat( Edit_dt->Text );         // 時間刻み
	total = StrToInt( Edit_number->Text );    // 地震波個数
	if ( CheckBox_win->Checked )
		band  = StrToFloat( Edit_haba->Text );    // バンド幅
	else
		band  = 0. ;

	accr = new double[total];
	i = 0;
	while( fscanf(fpin,"%f",&z ) != EOF ){
		accr[i] = (double) z;
		i++ ;
		if ( i >= total ) break;
	}
	fclose(fpin);
	n = i;
	a = log( (double)(n-1) ) / log( 2. ) + 1.0 ;
	nn = pow ( 2., (int) a ) ;
	nfold = nn / 2 ;
	df = 1. / nn / dt ;

	time = new double[nn];
	freq = new double[nn];
	afft = new double[nn];
	comp = new struct complex<double>[nn];
	pow_spec = new double[nn];
	temp1    = new double[nn];
	temp2    = new double[nn];

	for ( i=0; i<nn; i++ ) {
		time[i] = dt * i ;
		freq[i] = (double)i / nn / dt ;
	}
                              // FFTの準備
	for ( i=0; i<n; i++ )
		comp[i] = complex<double>( accr[i], 0. ) ;
	for ( i=n; i<nn; i++ )
		comp[i] = zero ;
                              // フーリエ変換で加速度のフーリエスペクトル
	fast  ( comp, nn, -1 ) ;
	for ( i=0; i<=nfold; i++ ) {
		afft[i] = abs( comp[i] ) * dt ;
	}

                              // フーリエスペクトルの平滑化
	if ( smooth( nfold, afft, pow_spec, df, band, temp1, temp2 ) < 0 )
		return ;

// ファイルに出力

	fprintf(fpout,"%s\t%s\t%s\n",
		"周波数(Hz)", "フーリエスペクトル", "パワースペクトル" );
	for ( i=0; i<=nfold; i++ )
		fprintf(fpout,"%lf\t%lf\t%lf\n",
			freq[i], afft[i], pow_spec[i] );

	fclose(fpout);

	delete[] time ;
	delete[] freq ;
	delete[] accr ;
	delete[] temp1 ;
	delete[] temp2 ;
	delete[] afft ;
	delete[] pow_spec ;

	Application->MessageBox("計算が終了しました。",
		"スペクトル解析",MB_ICONEXCLAMATION | MB_OK ) ;
}