package DSP::LinPred;
use 5.008005;
use Mouse;
our $VERSION = "0.05";

has 'mu' => (
    is => 'rw',
    isa => 'Num',
    default => 0.001
    );
has 'mu_mode' => (
    is => 'rw',
    isa => 'Int',
    default => 0
    );
has 'h_length' => (
    is => 'rw',
    isa => 'Int',
    default => 100
    );
has 'h' => (
    is => 'rw',
    isa => 'ArrayRef[Num]',
    default => sub{[(0) x 100]}
    );
has 'x_stack' => (
    is => 'rw',
    isa => 'ArrayRef[Num]',
    default => sub{[(0) x 100]}
    );
has 'x_count' => (
    is => 'rw',
    isa => 'Int',
    default => 0
    );
has 'current_error' => (
    is => 'rw',
    isa => 'Num',
    default => 0
    );
has 'dc' => (
    is => 'rw',
    isa => 'Num',
    default => 0
    );
has 'dc_mode' => (
    is => 'rw',
    isa => 'Int',
    default => 1
    );
has 'dc_init' => (
    is => 'rw',
    isa => 'Num',
    default => 0
    );
has 'stddev' => (
    is => 'rw',
    isa => 'Num',
    default => 0
    );
has 'stddev_mode' => (
    is => 'rw',
    isa => 'Int',
    default => 1
    );
has 'stddev_init' => (
    is => 'rw',
    isa => 'Num',
    default => 1
    );
has 'iir_mode' => (
    is => 'rw',
    isa => 'Int',
    default => 0
    );
has 'iir_a' => (
    is => 'rw',
    isa => 'Num',
    default => 0.01
    );



# filter specification
# mu : step size
# h_length : filter size
sub set_filter{
    my $self = shift;
    my $conf = shift;
    if(defined($conf->{mu_mode})){
        $self->mu_mode($conf->{mu_mode});
    }
    if(defined($conf->{iir_mode})){
        $self->iir_mode($conf->{iir_mode});
    }
    if(defined($conf->{iir_a})){
        $self->iir_a($conf->{iir_a});
    }
    if(defined($conf->{filter_length})){
	$self->h_length($conf->{filter_length});
	$self->h([(0) x $conf->{filter_length}]);
	if(defined($conf->{dc_init})){
	    $self->x_stack([($conf->{dc_init}) x $conf->{filter_length}]);
	}else{
	    $self->x_stack([(0) x $conf->{filter_length}]);
	}
	if(defined($conf->{iir_a})){
	    $self->iir_a($conf->{iir_a});
	}else{
	    $self->iir_a(1 / $conf->{filter_length})
	}
    }
    if(defined($conf->{dc_mode})){
	$self->dc_mode($conf->{dc_mode});
    }
    if(defined($conf->{dc_init})){
	$self->dc($conf->{dc_init});
	$self->dc_init($conf->{dc_init});
    }
    if(defined($conf->{stddev_mode})){
        $self->stddev_mode($conf->{stddev_mode});
    }
    if(defined($conf->{stddev_init})){
        $self->stddev($conf->{stddev_init});
        $self->stddev_init($conf->{stddev_init});
    }
}

# reset filter state
sub reset_state{
    my $self = shift;
    my $h_length = $self->h_length;
    $self->h([(0) x $h_length]);
    $self->x_stack([($self->dc_init) x $h_length]);
    $self->current_error(0);
    $self->dc($self->dc_init);
    $self->x_count(0);
    $self->stddev($self->stddev_init);
}

# prediction only
# predict_num : number of output predicted values
# this method returns list reference of predicted values
sub predict{
    my $self = shift;
    my $predict_num = shift;
    my $h = $self->h;
    my $x_stack = $self->x_stack;
    my $estimated;
    for(1 .. $predict_num){
	my $x_est = 0;
	for( my $k = 0; $k <= $#{$h} and $k <= $self->x_count; $k++){
	    $x_est += $h->[$k] * ($x_stack->[$k] - $self->dc);
	}
	$x_est += $self->dc;
	unshift(@$x_stack,$x_est);
	push(@$estimated,$x_est);
	pop(@$x_stack);
    }
    return($estimated);
}

# update only
# x should be array reference

sub update{
    my $self = shift;
    my $x = shift;
    my $h_length = $self->h_length;
    my $h = $self->h;
    my $x_stack = $self->x_stack;

    for ( my $kx=0; $kx <= $#{$x}; $kx++){
	unshift(@$x_stack,$x->[$kx]);
	pop(@$x_stack);
	$self->x_count($self->x_count + 1);
	if($self->dc_mode == 1){
	    if($self->iir_mode == 1){
		$self->dc_update_iir($x->[$kx]);
	    }else{
		$self->dc_update;
	    }
	}
	if($self->stddev_mode == 1){
	    if($self->iir_mode == 1){
		$self->stddev_update_iir($x->[$kx]);
	    }else{
		$self->stddev_update;
	    }
	}
	my $x_est = 0;
	for( my $k = 0; $k <= $#{$h} and $k <= $self->x_count;$k++){
	    $x_est += $h->[$k] * ($x_stack->[$k] - $self->dc);
	}
	my $error = $x->[$kx] - ($x_est + $self->dc);
	$self->current_error($error);
	my $h_new = $h;
	my $tmp_coef = 1;
	if($self->stddev_mode == 1){
	    $tmp_coef = $self->mu * $error / (1 + $self->stddev);
	}else{
	    $tmp_coef = $self->mu * $error;
	}
	if($self->mu_mode == 1){
	    $tmp_coef = 10 * $self->mu / (1 + $self->h_length);
	}

	for(my $k = 0;$k <= $#{$h} and $k <= $self->x_count; $k++){
	    $h_new->[$k] = 
		$h->[$k] 
		+ $tmp_coef * ($x_stack->[$k] - $self->dc);
	}
	$self->h($h_new);
    }
}

## DC component calculation and update
# using x_stack

sub dc_update{
    my $self = shift;
    my $x_stack = $self->x_stack;
    my $mean = 0;
    my $num = $#$x_stack + 1;
    for(0 .. $#$x_stack){
        $mean += $x_stack->[$_];
    }
    $mean = $mean / $num;
    $self->dc($mean);
}

## update by IIR
sub dc_update_iir{
    my $self = shift;
    my $x = shift;
    my $new_dc = ($x - $self->dc * $self->iir_a)/(1 - $self->iir_a);
    $self->dc($new_dc);
}


## standard deviation calculation and update
sub stddev_update{
    my $self = shift;
    my $x_stack = $self->x_stack;
    my $variance = 0;
    my $num = $#$x_stack + 1;
    for(0 .. $#$x_stack){
        $variance += ($x_stack->[$_] - $self->dc)**2;
    }
    $variance = $variance / $num;
    $self->stddev(sqrt($variance));
}

## update by IIR
sub stddev_update_iir{
    my $self = shift;
    my $x = shift;
    my $tmp = sqrt(($x - $self->dc)**2);
    my $new_stddev = ($tmp - $self->stddev * $self->iir_a)/(1 - $self->iir_a);
    $self->stddev($new_stddev);
}

## calculation of mean value of filter
sub filter_dc{
    my $self = shift;
    my $h = $self->h;
    my $mean = 0;
    my $num = $#$h + 1;
    for(0 .. $#$h){
        $mean += $h->[$_];
    }
    return($mean / $num);
}

## calculation of stddev of filter
sub filter_stddev{
    my $self = shift;
    my $h = $self->h;
    my $variance = 0;
    my $num = $#$h + 1;
    for(0 .. $#$h){
        $variance += ($h->[$_])**2;
    }
    return(sqrt($variance / $num));
}



1;
__END__

=encoding utf-8

=head1 NAME

DSP::LinPred - Linear Prediction

=head1 SYNOPSIS

    use DSP::LinPred;

    # OPTIONS
    # mu       : Step size of filter. (default = 0.001)
    #
    # h_length : Filter size. (default = 100)
    # dc_mode  : Direct Current Component estimation.
    #            it challenges to estimating DC component when set 1.
    #            (default = 1 enable)
    # dc_init  : Initial DC bias.
    #            It *SHOULD* be set value *ACCURATELY* when dc_mode => 0.
    #            (default = 0)
    #
    # stddev_mode : Step size correction by stddev of input.
    #               (default = 1 enable)
    # stddev_init : Initial value of stddev.
    #               (default = 1)
    #

    my $lp = DSP::LinPred->new;

    # set filter
    $lp->set_filter({
                     mu => 0.001,
                     filter_length => 500,
                     dc_mode => 1,
                     stddev_mode => 1
                    });

    # defining signal x
    my $x = [0,0.1,0.5, ... ]; # input signal

    # Updating Filter
    $lp->update($x);
    my $current_error = $lp->current_error; # get error

    # Prediction
    my $pred_length = 10;
    my $pred = $lp->predict($pred_length);
    for( 0 .. $#$pred ){ print $pred->[$_], "\n"; }


=head1 DESCRIPTION

DSP::LinPred is Linear Prediction by Least Mean Squared Algorithm.

This Linear Predicting method can estimate the standard deviation, direct current component, and predict future value of input.

=head1 METHODS

=head2 I<set_filter>

I<set_filter> method sets filter specifications to DSP::LinPred object.

    $lp->set_filter(
        {
            mu => $step_size, # <Num>
            filter_length => $filter_length, # <Int>
            dc_init => $initial_dc_bias, # <Num>
            dc_mode => $dc_estimation, # <Int>, enable when 1
            stddev_init => $initial_stddev, # <Num>
            stddev_mode => $stddev_estimation # <Int>, enable when 1
        });


=head2 I<update>

I<update> method updates filter state by source inputs are typed ArrayRef[Num].

    my $x = [0.13,0.3,-0.2,0.5,-0.07];
    $lp->update($x);

If you would like to extract the filter state, you can access member variable directly like below.

    my $filter = $lp->h;
    for( 0 .. $#$filter ){ print $filter->[$_], "\n"; }

=head2 I<predict>

I<predict> method generates predicted future values of inputs by filter.

    my $predicted = $lp->predict(7);
    for( 0 .. $#$predicted ){ print $predicted->[$_], "\n";}

=head2 I<filter_dc>

This method can calculate mean value of current filter.

    my $filter_dc = $lp->filter_dc;

=head2 I<filter_stddev>

This method can calculate standard deviation of current filter.

    my $filter_stddev = $lp->filter_stddev;


=head1 READING STATES

=head2 I<current_error>

    # It returns value of current prediction error
    # error = Actual - Predicted
    my $current_error = $lp->current_error;
    print 'Current Error : '.$current_error, "\n";

=head2 I<h>

    # It returns filter state(ArrayRef)
    my $filter = $lp->h;
    print "Filter state\n";
    for( 0 .. $#$filter ){ print $_.' : '.$filter->[$_],"\n"; }

=head2 I<x_count>

    # It returns value of input counter used in filter updating.
    my $x_count = $lp->x_count;
    print 'Input count : '.$x_count, "\n";

=head2 I<dc>

    # Get value of current Direct Current Components of inputs.
    my $dc = $lp->dc;
    print 'Current DC-Component : '.$dc, "\n";

=head2 I<stddev>

    # Get value of current standard deviation of inputs.
    my $stddev = $lp->dc;
    print 'Current STDDEV : '.$stddev, "\n";

=head1 LICENSE

Copyright (C) Toshiaki Yokoda.

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself.

=head1 AUTHOR

Toshiaki Yokoda E<lt>E<gt>

=cut