1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
|
=head1 NAME
PDL::Filter::Linear - linear filtering for PDL
=head1 SYNOPSIS
$a = new PDL::Filter::Linear(
{Weights => $v,
Point => 10});
$b = new PDL::Filter::Gaussian(15,2); # 15 points, 2 std devn.
($pred,$corrslic) = $a->predict($dat);
=head1 DESCRIPTION
A wrapper for generic linear filters.
Just for convenience. This should in the future use DataPresenter.
Also, this class should at some point learn to do FFT whenever it is
useful.
=cut
package PDL::Filter::Linear;
use PDL;
use PDL::Basic;
use PDL::Slices;
use PDL::Primitive;
use strict;
sub new($$) {
my($type,$pars) = @_;
my $this = bless {},$type;
$this->{Weights} = ((delete $pars->{Weights}) or
barf("Must specify weights\n"));
$this->{Point} = (delete $pars->{Point} or 0);
$this;
}
sub predict($$) {
my($this,$data) = @_;
my $ldata = $data->lags(0,1,$this->{Weights}->getdim(0));
inner($ldata->xchg(0,1),$this->{Weights},
(my $pred = PDL->null));
return wantarray ? ($pred,$ldata->slice(":,($this->{Point})")) :
$pred ;
}
package PDL::Filter::Gaussian;
use PDL; use PDL::Basic; use PDL::Slices; use PDL::Primitive;
use strict;
@PDL::Filter::Gaussian::ISA = qw/PDL::Filter::Linear/;
sub new($$) {
my($type,$npoints,$sigma) = @_;
my $cent = int($npoints/2);
my $x = ((PDL->zeroes($npoints )->xvals) - $cent)->float;
my $y = exp(-($x**2)/(2*$sigma**2));
# Normalize to unit total
$y /= sum($y);
return PDL::Filter::Linear::new($type,{Weights => $y,
Point => $cent});
}
# Savitzky-Golay (see Numerical Recipes)
package PDL::Filter::SavGol;
use PDL; use PDL::Basic; use PDL::Slices; use PDL::Primitive;
use strict;
@PDL::Filter::Gaussian::ISA = qw/PDL::Filter::Linear/;
# XXX Doesn't work
sub new($$) {
my($type,$deg,$nleft,$nright) = @_;
my $npoints = $nright + $nleft + 1;
my $x = ((PDL->zeroes($npoints )->xvals) - $nleft)->float;
my $mat1 = ((PDL->zeroes($npoints,$deg+1)->xvals))->float;
for(0..$deg-1) {
(my $tmp = $mat1->slice(":,($_)")) .= ($x ** $_);
}
my $y;
# Normalize to unit total
return PDL::Filter::Linear::new($type,{Weights => $y,
Point => $nleft});
}
|
