package MEM;
use Exporter;
@ISA = qw(Exporter);

#公開したいサブルーチン
@EXPORT = qw();

use strict;

use Sci qw($pi tanh cosh sinh);

my $pi          = 3.1415926535;

#================================================
# 静的メンバー関数
#================================================
sub MEM
{
	my ($dt, $px, $mmx, $lmax) = @_;

	my $nmax = scalar @$px;

	my ($i , $lag , $lg1 , $lx1);
	my (@f, @fpe, @spec, @c, @g);

	my @x2 = (0.0, @$px);
	my $pm = MEM::MEMBurg($nmax , $mmx , $lmax , \@x2 , \@g , \@c , \@fpe);
#	my $pm = MEM::MEMBurg($nmax , $mmx , $lmax , $px , \@g , \@c , \@fpe);
	MEM::MEMSub($nmax , $mmx , $dt , \@g , \@f , \@spec , $pm);

	return (\@f, \@fpe, \@spec, \@c, $mmx, $lmax, $pm);
}

sub MEMSub
{
	my ($nmax , $mmx , $dt , $pg , $pf , $pspec , $pm) = @_;

	my $f0 = 1.0 / $dt / ($nmax - 1.0);
	for(my $i = 1 ; $i <= $nmax ; $i++)
#	for(my $i = 0 ; $i < $nmax ; $i++)
	{
		$pf->[$i] = $f0 * $i;
		my $sumr = 1.0;
		my $sumi = 0.0;
		my $cr = 0.0;
		my $ci = 2.0 * $pi * $pf->[$i] * $dt;
		for(my $j = 1 ; $j <= $mmx ; $j++)
#		for(my $j = 0 ; $j < $mmx ; $j++)
		{
			$sumr += $pg->[$j] * cos($j * $ci);
			$sumi += $pg->[$j] * sin($j * $ci);
		}
		$pspec->[$i] = 2.0 * $pm * $dt / ($sumr * $sumr + $sumi * $sumi);
	}
}

sub MEMBurg
{
	my ($nmax , $mmx , $lmax , $px , $pg , $pc , $pfpe) = @_;

	my (@b1, @b2, @a);

	my $sum = 0.0;
	for(my $i = 1 ; $i <= $nmax ; $i++) {
#	for(my $i = 0 ; $i < $nmax ; $i++) {
		$sum += $px->[$i] * $px->[$i];
	}

	$pc->[1] = $sum / $nmax;
#	$pc->[0] = $sum / $nmax;
	my $pm = $pc->[1];
#	my $pm = $pc->[0];
#	$pfpe->[0] = ($nmax + 1.0) / ($nmax - 1.0) * $pm;
	$pfpe->[1] = ($nmax + 1.0) / ($nmax - 1.0) * $pm;

	my $lag = $mmx + 1;
	my $lg1 = $lag + 1;
	$b1[1] = $px->[1];
#	$b1[0] = $px->[0];
	for(my $i = 2 ; $i <= $nmax ; $i++)
#	for(my $i = 1 ; $i < $nmax ; $i++)
	{
		$b1[$i]   = $px->[$i];
		$b2[$i-1] = $px->[$i];
	}
	
	for(my $m = 1 ; $m <= $mmx ; $m++)
#	for(my $m = 0 ; $m < $mmx ; $m++)
	{
#printf("\n M=%4d " , $m);
		MEM::MEMLevins($nmax , $pg , $pm , \@b1 , \@b2 , \@a , $m);
		$sum = 0.0;
		for(my $i = 1 ; $i <= $m ; $i++) {
#		for(my $i = 0 ; $i < $m ; $i++) {
			$sum -= $pc->[$m-$i] * $pg->[$i];
		}
		$pc->[$m+1] = $sum;
		if($m != $nmax-1) {
			$pfpe->[$m+1] = ($nmax + $m + 1) / ($nmax - $m - 1) * $pm;
		}
	}
	
	if($lag < $lmax) {
		for(my $l = $lg1 ; $l <= $lmax ; $l++)
		{
			$sum = 0.0;
			for(my $i = 1 ; $i <= $mmx ; $i++) {
#			for(my $i = 0 ; $i < $mmx ; $i++) {
				$sum -= $pc->[$l] * $pg->[$i];
			}
			$pc->[$l] = $sum;
		}
	}

	return $pm;
}

sub MEMLevins
{
	my ($nmax , $pg , $pm , $pb1 , $pb2 , $pa , $m) = @_;

	my ($sn, $sd) = (0.0, 0.0);
	for(my $i = 1 ; $i <= $nmax - $m ; $i++)
#	for(my $i = 0 ; $i < $nmax - $m ; $i++)
	{
		$sn += $pb1->[$i] * $pb2->[$i];
		$sd += $pb1->[$i] * $pb1->[$i] + $pb2->[$i] * $pb2->[$i];
	}
	$pg->[$m] = -2.0 * $sn / $sd;
	$pm *= (1.0 - $pg->[$m] * $pg->[$m]);
	if($m != 1) {
		for(my $k = 1 ; $k <= $m-1 ; $k++) {
#		for(my $k = 0 ; $k < $m-1 ; $k++) {
			$pg->[$k] = $pa->[$k] + $pg->[$m] * $pa->[$m-$k];
		}
	}
	for(my $i = 1 ; $i <= $m-1 ; $i++) {
#	for(my $i = 0 ; $i < $m ; $i++) {
		$pa->[$i] = $pg->[$i];
	}
	for(my $i = 1 ; $i <= $nmax - $m ; $i++)
#	for(my $i = 0 ; $i < $nmax - $m - 1 ; $i++)
	{
		$pb1->[$i] += $pa->[$m] * $pb2->[$i];
		$pb2->[$i] = $pb2->[$i+1] + $pa->[$m] * $pb1->[$i+1];
	}
}

sub MEM2
{
	my ($dt, $px, $isw, $mmax) = @_;
# $dt  : サンプル間隔
# $isw : 1: Model order is given, 0: Model order is determined with the upper limitation of $mmax
# $mmax: モデル次数 あるいは 最大モデル次数
#printf(" ***  ｻﾝﾌﾟﾙ ﾃﾝｽｳ       =");
	my $nd    = scalar @$px;
	my $ns    = $nd;
	my $wnmin = 0.0;
	my $wnmax = $nd;
	my $ci    = 2.0 * $pi;

	my (@fpe, @rfpe, @r, @rr);
	my ($minm, $pmm) = MEM::gs1420($nd, $px, $isw, $mmax, \@fpe, \@rfpe, \@r, \@rr); # Burg法による自己回帰係数と予測誤差の推定

	my $wnmin = 0.0;
	my $wnmax = 1.0 / $dt;
# スペクトル点数
	my $ns = $nd;
	my $wnint = ($wnmax - $wnmin) / $ns; # スペクトルの計算
	MEM::gs1960($dt, $px, $ci, $minm, $pmm, $wnmin, $wnmax, $wnint, \@rfpe); # スペクトルの計算

	my @f;
	for(my $i = 1 ; $i <= $ns ; $i++) {
#	for(my $i = 0 ; $i < $ns ; $i++) {
		$f[$i] = $wnint * $i;
	}

	return (\@f, \@fpe, \@rfpe, \@r, $mmax, $mmax, $pmm);
}

# スペクトルの計算
sub gs1960 {
	my ($dt, $px, $ci, $minm, $pmm, $wnmin, $wnmax, $wnint, $prfpe) = @_;

	$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 += $prfpe->[$j] * cos($ci * $f * $j);
			$sum2 += $prfpe->[$j] * sin($ci * $f * $j);
		}
		$px->[$i] = $pmm * $dt / ($sum1 * $sum1 + $sum2 * $sum2);
		$i = $i + 1;
	}
}


# Burg法による自己回帰係数と予測誤差の推定
sub gs1420 {
	my ($nd, $px, $isw, $mmax, $pfpe, $prfpe, $pr, $prr) = @_;

	my ($minm, $pmm);
	my ($z, $m);

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

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

	my $fpemin = ($nd + 1.0) / ($nd - 1.0) * $pm;
	$pfpe->[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 += $px->[$i] * $y[$i];
			$sumd += $px->[$i] * $px->[$i] + $y[$i] * $y[$i];
		}
		my $rm = -2.0 * $sumn / $sumd;
		$pr->[$m] = $rm;
		$pm *= 1.0 - $rm * $rm;
		if($m != 1) {
			for(my $i = 0 ; $i < $m - 1 ; $i++)
			{
				$pr->[$i] = $prr->[$i] + $rm * $prr->[$m-$i];
			}
		}
		for(my $i = 0 ; $i < $m ; $i++)
		{
			$prr->[$i] = $pr->[$i];
		}
		for(my $i = 0 ; $i < $nd - $m - 1 ; $i++)
		{
			$px->[$i] = $px->[$i] + $rm * $y[$i];
			$y[$i]    = $y[$i+1]  + $rm * $px->[$i+1];
		}
		if($isw == 0) {
			($minm, $pmm) = MEM::gs1800($nd, $fpemin, $m, $pm, $pfpe, $prfpe, $pr);
		}
	}
	if($isw == 1) {
		($minm, $pmm) = MEM::gs1880($m, $pm, $pr, $prfpe);
	}

	return ($minm, $pmm);
}

sub gs1880
{
	my ($mmax, $pm, $pr, $prfpe) = @_;

	my $minm = $mmax;
	my $pmm = $pm;
	for(my $i = 0 ; $i < $mmax ; $i++)
	{
		$prfpe->[$i] = $pr->[$i];
	}

	return ($minm, $pmm);
}

sub gs1800
{
	my ($nd, $fpemin, $m, $pm, $pfpe, $prfpe, $pr) = @_;

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

1;