#!/usr/bin/perl

use strict;

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

use Utils;
use JFile;
use Sci::MEM;
use Sci qw($pi tanh cosh sinh);

my ($x0, $x1, $nData) = (0, 1.0, 128);
#my ($x0, $x1, $nData) = (0, 1.0, 1024);

my $xstep = ($x1 - $x0) / $nData;
print "x = ($x0, $x1, $xstep)\n";
print "nData=$nData\n";

my (@x, @y);
for(my $i = 0 ; $i < $nData ; $i++) {
	$x[$i] = $x0 + $i * $xstep;
	$y[$i] = &Func($x[$i]);
}
#my $df = 1.0 / ($pX->[1] + $pX->[$n-1] - 2.0 * $pX->[0]);

my $mmx  = 38;
my $lmax = 50;

&MEMCall(\@x, \@y, $mmx, $lmax);
exit;


sub MEMCall
{
	my ($pX, $pY, $mmx, $lmax) = @_;

	my $lag = $mmx + 1;
	if($lag+1 > $lmax) {
		$lmax = $lag + 1;
	}

	my $nmax = scalar @$pX;

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

	my ($lg1 , $lx1, $pm);

	my $dt = $x[1] - $x[0];
#	my ($pf, $pfpe, $pspec, $pc, $mmx, $lmax, $pm) = MEM::MEM($dt, $pY, $mmx, $lmax);
	my $isw  = 0;
	my $mmax = $mmx;
	my ($pf, $pfpe, $pspec, $pc, $mmx, $lmax, $pm) = MEM::MEM2($dt, $pY, $isw, $mmax);

	$lag = $mmx + 1;
	$lg1 = $lag + 1;
	$lx1 = $lmax + 1;
	printf("\nf\tspc\tfpe\tauto c.");
	for(my $i = 0 ; $i < $lag ; $i++) {
		printf("\n%g\t%7.3f\t%g\t%g" ,
				$pf->[$i] , $pspec->[$i] , $pfpe->[$i] , $pc->[$i]);
	}
	if($lag != $lmax) {
		for(my $i = $lg1 ; $i <= $lmax ; $i++) {
			printf("\n%g\t%g\t%g" , $pf->[$i] , $pspec->[$i] , $pc->[$i]);
		}
	}
	for(my $i = $lx1 ; $i < $nmax ; $i++) {
		printf("\n%g\t%g" , $pf->[$i] , $pspec->[$i]);
	}
}

sub Func
{
	my ($x) = @_;

$x += rand(0.03);
	
	my ($f1, $p1, $A1) = ( 1.5, $pi/4.0, 1.0);
	my ($f2, $p2, $A2) = ( 3.0, $pi/3.0, 0.3);
	my ($f3, $p3, $A3) = (10.0, $pi/6.0, 0.5);
	return $A1 * sin(2.0*$pi * $f1 * $x + $p1)
		 + $A2 * sin(2.0*$pi * $f2 * $x + $p2)
		 + $A3 * sin(2.0*$pi * $f3 * $x + $p3);
}