#!/usr/bin/perl

use strict;

use lib 'd:/Programs/Perl/lib';

use Utils;
use JFile;

my $MMAX = 300;
my $ci = 2*3.14159;

my (@x, @y, @fpe, @r, @rr, @rfpe);
my ($isw, $nd, $dt, $an, $mmax, $ns);
my ($wnmin, $wnmax, $wnint);
my ($fpemin);
my ($pm, $minm, $pmm);
#double rm , wnmin , wnmax , sum1 , sum2 , sum , sumn , sumd ,wnint;
#int    isw , nd , mmax , isw , ns , minm;
#double dt , pmm , pm , av;

&MEM();
exit;


sub MEM
{
	my ($an, @p, @n);

	printf("      ***** MEM *****\n\n");
	printf(" ***  ｻﾝﾌﾟﾙ ﾃﾝｽｳ       =");
	$nd = <> + 0;
	printf("nd:%d\n" ,$nd);
	printf(" ***   ｻﾝﾌﾟﾙ ｶﾝｶｸ       =");
	$dt = <> + 0.0;
	printf("dt:%g\n" , $dt);

	printf(" ﾓﾃﾞﾙ ｼﾞｽｳ ｦ ｱﾀｴﾏｽｶ (Y/N)\n");
	$an = <>;
	Utils::DelSpace($an);

	my $p;
	if($an eq 'Y' or $an eq 'y')
	{
		$isw = 1;
		$p = " ***   ﾓﾃﾞﾙ ｼﾞｽｳ           =";
	}
	elsif($an eq 'N' or $an eq 'n')
	{
		$isw = 0;
		$p = " ***   ｻｲﾀﾞｲ ﾓﾃﾞﾙ ｼﾞｽｳ     =";
	}
	else {
		print "Invalid answer [$an].\n";
		exit;
	}

	printf("%s", $p);
	$mmax = <> + 0;
	printf("mmax:%d\n", $mmax);

	&gs1340(); #データ入力
	$ns = $nd;
	$wnmin = 0.0;
	$wnmax = $nd;
	my $n = " Data ";
#/*	gs2110();*/ # データのグラフ表示
	&gs1420(); # Burg法による自己回帰係数と予測誤差の推定

	printf("   ｻｲｼｮｳ ｼｭｳﾊｽｳ =");
	$wnmin = <> + 0.0;
	printf("   ｻｲﾀﾞｲ ｼｭｳﾊｽｳ =");
	$wnmax = <> + 0.0;
	printf("wnmin,max:%g,  %g\n", $wnmin, $wnmax);
	printf("   ｽﾍﾟｸﾄﾙ ﾃﾝｽｳ  =");
	$ns = <> + 0;
	$wnint = ($wnmax - $wnmin) / $ns; # スペクトルの計算
	printf("ns:%d   wnint:%g\n" , $ns , $wnint);
	&gs1960(); # スペクトルの計算
	$n = " MEM Spectrum ";
	&gs2110(); # スペクトル推定結果のグラフ

	return 1;
}

# スペクトル推定結果のグラフ
sub gs2110
{
	my $out = JFile->new('out.txt', 'w') or die "$!: Can not write to [out.txt].\n";
	
	for(my $i = 0 ; $i < $ns ; $i++)
	{
		$out->printf("%g\n", $x[$i]);
	}
	$out->Close();
}

# スペクトルの計算
sub gs1960 {
	double f;

	$ci *= $dt;
	my $i = 0;
	for(my $f = $wnmin ; $f <= $wnmax ; $f += $wnint)
	{
		my $sum1 = 1.0;
		my $sum2 = 0.0;
		for(my $j = 0 ; $j < $minm ; $j++)
		{
			$sum1 += $rfpe[$j] * cos($ci * $f * $j);
			$sum2 += $rfpe[$j] * sin($ci * $f * $j);
		}
		$x[$i] = $pmm * $dt / ($sum1 * $sum1 + $sum2 * $sum2);
		$i = $i + 1;
	}
}

# Burg法による自己回帰係数と予測誤差の推定
sub gs1420 {
	my ($z, $m);

	my $sum = 0.0;
	for(my $i = 0 ; $i < $nd ; $i++) {
		$sum += $x[$i];
	}
	my $av = $sum / $nd;

	$sum = 0.0;
	for(my $i = 0; $i < $nd; $i++)
	{
		$z = $x[$i] - $av;
		$x[$i]   = $z;
		$y[$i-1] = $z;
		$sum += $z*$z;
	}
	my $pm = $sum / $nd;

	$fpemin = ($nd + 1.0) / ($nd - 1.0) * $pm;
	$fpe[0] = $fpemin;
	for(my $m = 0; $m < $mmax; $m++)
	{
		my $sumn = 0.0;
		my $sumd = 0.0;
		for(my $i = 0; $i < $nd - $m ; $i++)
		{
			$sumn += $x[$i] * $y[$i];
			$sumd += $x[$i] * $x[$i] + $y[$i] * $y[$i];
		}
		my $rm = -2.0 * $sumn / $sumd;
		$r[$m] = $rm;
		$pm *= 1.0 - $rm * $rm;
		if($m != 1) {
			for(my $i = 0 ; $i < $m - 1 ; $i++)
			{
				$r[$i] = $rr[$i] + $rm * $rr[$m-$i];
			}
		}
		for(my $i = 0 ; $i < $m ; $i++)
		{
			$rr[$i] = $r[$i];
		}
		for(my $i = 0 ; $i < $nd - $m - 1 ; $i++)
		{
			$x[$i] = $x[$i] + $rm * $y[$i];
			$y[$i] = $y[$i+1] + $rm * $x[$i+1];
		}
		if($isw == 0) {
			&gs1800($m);
		}
	}
	if($isw == 1) {
		&gs1880();
	}
}

sub gs1880
{
	$minm = $mmax;
	$pmm = $pm;
	for(my $i = 0 ; $i < $mmax ; $i++)
	{
		$rfpe[$i] = $r[$i];
	}
}

sub gs1800
{
	my ($m) = @_;

	$fpe[$m] = ($nd + $m + 1.0) / ($nd - $m - 1.0) * $pm;
	if($fpe[$m] > $fpemin) {
		return;
	}
	$fpemin = $fpe[$m];
	$minm = $m;
	$pmm = $pm;
	for(my $i = 0 ; $i < $m ; $i++)
	{
		$rfpe[$i] = $r[$i];
	}
}



#データ入力
sub gs1340
{
	my @nm;
	my $fp;
	my $i;
	my ($ccc, $bbb);
	
	printf(" ***   ﾃﾞｰﾀﾌｧｲﾙﾒｲ       =");
	my $nm = <>;
	Utils::DelSpace($nm);
	my $in = JFile->new($nm, 'r') or die "$!: Can not read [$nm].\n";

	for(my $i = 0; my $i < $nd; $i++)
	{
		my $line = $in->ReadLine();
		($x[$i], $bbb, $ccc) = Utils::Split("\\s+", $line);
		printf("I=%d  X=%g\n" , $i , $x[$i]);
	}
	$in->Close();
}