#!/usr/bin/perl

use strict;

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

use Utils;
use JFile;

my $PI = 3.14159265358979323846;
my (@x, @f, @fpe, @spec, @g, @c);

&MEM();
exit;


sub MEM
{
	my ($nmax , $mmx , $lmax , $i , $lag , $lg1 , $lx1);
	my ($dt , $sum , $av , $pm);
	
	($nmax , $mmx , $lmax , $dt) = &mk_dat(\@x , $nmax , $mmx , $lmax , $dt);

	$lag = $mmx + 1;
	if($lag+1 > $lmax) {
		$lmax = $lag + 1;
	}
	
	$sum = 0.0;
	for($i = 0 ; $i < $nmax ; $i++) {
		$sum += $x[$i];
	}
	$av = $sum / $nmax;
	for($i = 0 ; $i < $nmax ; $i++) {
		$x[$i] -= $av;
	}
	
	$pm = &burg($nmax , $mmx , $lmax , \@x , \@g , \@c , \@fpe);
	&mem_1($nmax , $mmx , $dt , \@g , \@f , \@spec , $pm);
	
	$lag = $mmx + 1;
	$lg1 = $lag + 1;
	$lx1 = $lmax + 1;
	printf("\nfrec\tspc\tfpe\tauto c.");
	for($i = 0 ; $i < $lag ; $i++) {
		printf("\n%6.1f\t%7.3f\t%7.3f\t%7.3f" ,
				$f[$i] , $spec[$i] , $fpe[$i] , $c[$i]);
	}
	if($lag != $lmax) {
		for($i = $lg1 ; $i <= $lmax ; $i++) {
			printf("\n%6.1f\t%7.3f\t%6.3f" , $f[$i] , $spec[$i] , $c[$i]);
		}
	}
	for($i = $lx1 ; $i < $nmax ; $i++) {
		printf("\n%6.1f\t%7.3f" , $f[$i] , $spec[$i]);
	}
}

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

	my $f0 = 1.0 / (($nmax - 1.0) * $dt);
	for(my $i = 0 ; $i < $nmax ; $i++)
	{
		$pf->[$i] = $f0 * $i;
		printf("\n  freq=%6.2f" , $pf->[$i]);
		my $sumr = 1.0;
		my $sumi = 0.0;
		my $cr = 0.0;
		my $ci = 2.0 * $PI * $pf->[$i] * $dt;
		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 burg
{
	my ($nmax , $mmx , $lmax , $px , $pg , $pc , $pfpe) = @_;

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

	my $sum = 0.0;
	for(my $i = 0 ; $i < $nmax ; $i++) {
		$sum += $px->[$i] * $px->[$i];
	}
	$c[0] = $sum / $nmax;
	my $pm = $c[0];
	$fpe[0] = ($nmax + 1.0) / ($nmax - 1.0) * $pm;
	my $lag = $mmx + 1;
	my $lg1 = $lag + 1;
	$b1[0] = $px->[0];
	for(my $i = 1 ; $i <= $nmax-1 ; $i++)
	{
		$b1[$i]   = $px->[$i];
		$b2[$i-1] = $px->[$i];
	}
	
	for(my $m = 0 ; $m < $mmx ; $m++)
	{
		printf("\n M=%4d " , $m);
		&levins($nmax , $pg , $pm , \@b1 , \@b2 , \@a , $m);
		$sum = 0.0;
		for(my $i = 0 ; $i < $m ; $i++) {
			$sum -= $pc->[$m-$i+1] * $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 = 0 ; $i < $mmx ; $i++) {
				$sum -= $pc->[$l-1] * $pg->[$i];
			}
			$pc->[$l] = $sum;
		}
	}

	return $pm;
}

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

	my ($sn, $sd) = (0.0, 0.0);
	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 = 0 ; $k < $m-1 ; $k++) {
			$pg->[$k] = $pa->[$k] + $pg->[$m] * $pa->[$m-$k];
		}
	}
	for(my $i = 0 ; $i < $m ; $i++) {
		$pa->[$i] = $pg->[$i];
	}
	for(my $i = 0 ; $i < $nmax - $m ; $i++)
	{
		$pb1->[$i] += $pa->[$m] * $pb2->[$i];
		$pb2->[$i] = $pb2->[$i+1] + $pa->[$m] * $pb1->[$i+1];
	}
}

sub mk_dat
{
	my ($px, $nmax, $mmx, $lmax, $dt) = @_;

	$nmax = 241;
	$mmx  = 38;
	$lmax = 50;
	$dt = 1.0 / ($nmax - 1);
	my $t1 = 1.0 / 124.5;
	my $t2 = 1.0 /  65.0;
	my $t3 = 1.0 /  64.0;
	my $t4 = 1.0 /   4.0;
	for(my $i = 0 ; $i < $nmax ; $i++)
	{
		my $t = $dt * $i;
		my $y = sin(2.0 * $PI * $t / $t1) + sin(2.0 * $PI * $t / $t2)
		      + sin(2.0 * $PI * $t / $t3) + sin(2.0 * $PI * $t / $t4);
		$px->[$i] = $y;
#		/* + 2.0 * drand48()  - 2;*/
	}
	return ($nmax , $mmx , $lmax , $dt);
}