File Coverage

blib/lib/DSP/LinPred.pm
Criterion Covered Total %
statement 97 130 74.6
branch 18 32 56.2
condition 8 9 88.8
subroutine 8 12 66.6
pod 5 10 50.0
total 136 193 70.4


line stmt bran cond sub pod time code
1             package DSP::LinPred;
2 2     2   31394 use 5.008005;
  2         8  
  2         77  
3 2     2   1584 use Mouse;
  2         67896  
  2         12  
4             our $VERSION = "0.06";
5              
6             has 'mu' => (
7             is => 'rw',
8             isa => 'Num',
9             default => 0.001
10             );
11             has 'mu_mode' => (
12             is => 'rw',
13             isa => 'Int',
14             default => 0
15             );
16             has 'h_length' => (
17             is => 'rw',
18             isa => 'Int',
19             default => 100
20             );
21             has 'h' => (
22             is => 'rw',
23             isa => 'ArrayRef[Num]',
24             default => sub{[(0) x 100]}
25             );
26             has 'x_stack' => (
27             is => 'rw',
28             isa => 'ArrayRef[Num]',
29             default => sub{[(0) x 100]}
30             );
31             has 'x_count' => (
32             is => 'rw',
33             isa => 'Int',
34             default => 0
35             );
36             has 'current_error' => (
37             is => 'rw',
38             isa => 'Num',
39             default => 0
40             );
41             has 'dc' => (
42             is => 'rw',
43             isa => 'Num',
44             default => 0
45             );
46             has 'dc_mode' => (
47             is => 'rw',
48             isa => 'Int',
49             default => 1
50             );
51             has 'dc_init' => (
52             is => 'rw',
53             isa => 'Num',
54             default => 0
55             );
56             has 'stddev' => (
57             is => 'rw',
58             isa => 'Num',
59             default => 0
60             );
61             has 'stddev_mode' => (
62             is => 'rw',
63             isa => 'Int',
64             default => 1
65             );
66             has 'stddev_init' => (
67             is => 'rw',
68             isa => 'Num',
69             default => 1
70             );
71             has 'iir_mode' => (
72             is => 'rw',
73             isa => 'Int',
74             default => 0
75             );
76             has 'iir_a' => (
77             is => 'rw',
78             isa => 'Num',
79             default => 0.01
80             );
81              
82              
83              
84             # filter specification
85             # mu : step size
86             # h_length : filter size
87             sub set_filter{
88 2     2 1 21 my $self = shift;
89 2         3 my $conf = shift;
90 2 50       8 if(defined($conf->{mu_mode})){
91 0         0 $self->mu_mode($conf->{mu_mode});
92             }
93 2 100       6 if(defined($conf->{iir_mode})){
94 1         6 $self->iir_mode($conf->{iir_mode});
95             }
96 2 50       9 if(defined($conf->{iir_a})){
97 0         0 $self->iir_a($conf->{iir_a});
98             }
99 2 50       8 if(defined($conf->{filter_length})){
100 2         25 $self->h_length($conf->{filter_length});
101 2         71 $self->h([(0) x $conf->{filter_length}]);
102 2 50       12 if(defined($conf->{dc_init})){
103 2         48 $self->x_stack([($conf->{dc_init}) x $conf->{filter_length}]);
104             }else{
105 0         0 $self->x_stack([(0) x $conf->{filter_length}]);
106             }
107 2 50       12 if(defined($conf->{iir_a})){
108 0         0 $self->iir_a($conf->{iir_a});
109             }else{
110 2         9 $self->iir_a(1 / $conf->{filter_length})
111             }
112             }
113 2 50       9 if(defined($conf->{dc_mode})){
114 2         6 $self->dc_mode($conf->{dc_mode});
115             }
116 2 50       8 if(defined($conf->{dc_init})){
117 2         7 $self->dc($conf->{dc_init});
118 2         7 $self->dc_init($conf->{dc_init});
119             }
120 2 50       7 if(defined($conf->{stddev_mode})){
121 2         6 $self->stddev_mode($conf->{stddev_mode});
122             }
123 2 50       8 if(defined($conf->{stddev_init})){
124 0         0 $self->stddev($conf->{stddev_init});
125 0         0 $self->stddev_init($conf->{stddev_init});
126             }
127             }
128              
129             # reset filter state
130             sub reset_state{
131 1     1 0 1506 my $self = shift;
132 1         6 my $h_length = $self->h_length;
133 1         36 $self->h([(0) x $h_length]);
134 1         28 $self->x_stack([($self->dc_init) x $h_length]);
135 1         14 $self->current_error(0);
136 1         4 $self->dc($self->dc_init);
137 1         3 $self->x_count(0);
138 1         6 $self->stddev($self->stddev_init);
139             }
140              
141             # prediction only
142             # predict_num : number of output predicted values
143             # this method returns list reference of predicted values
144             sub predict{
145 2     2 1 20 my $self = shift;
146 2         5 my $predict_num = shift;
147 2         7 my $h = $self->h;
148 2         8 my $x_stack = $self->x_stack;
149 2         4 my $estimated;
150 2         8 for(1 .. $predict_num){
151 2         5 my $x_est = 0;
152 2   66     5 for( my $k = 0; $k <= $#{$h} and $k <= $self->x_count; $k++){
  602         2830  
153 600         1391 $x_est += $h->[$k] * ($x_stack->[$k] - $self->dc);
154             }
155 2         7 $x_est += $self->dc;
156 2         3 unshift(@$x_stack,$x_est);
157 2         7 push(@$estimated,$x_est);
158 2         6 pop(@$x_stack);
159             }
160 2         9 return($estimated);
161             }
162              
163             # update only
164             # x should be array reference
165              
166             sub update{
167 2     2 1 3013 my $self = shift;
168 2         4 my $x = shift;
169 2         6 my $h_length = $self->h_length;
170 2         6 my $h = $self->h;
171 2         6 my $x_stack = $self->x_stack;
172              
173 2         4 for ( my $kx=0; $kx <= $#{$x}; $kx++){
  2002         10579  
174 2000         4520 unshift(@$x_stack,$x->[$kx]);
175 2000         2600 pop(@$x_stack);
176 2000         7884 $self->x_count($self->x_count + 1);
177 2000 100       6431 if($self->dc_mode == 1){
178 1000 50       2792 if($self->iir_mode == 1){
179 0         0 $self->dc_update_iir($x->[$kx]);
180             }else{
181 1000         2136 $self->dc_update;
182             }
183             }
184 2000 50       6741 if($self->stddev_mode == 1){
185 2000 50       5775 if($self->iir_mode == 1){
186 0         0 $self->stddev_update_iir($x->[$kx]);
187             }else{
188 2000         4451 $self->stddev_update;
189             }
190             }
191 2000         3027 my $x_est = 0;
192 2000   100     3361 for( my $k = 0; $k <= $#{$h} and $k <= $self->x_count;$k++){
  512898         2336129  
193 510898         1198281 $x_est += $h->[$k] * ($x_stack->[$k] - $self->dc);
194             }
195 2000         5196 my $error = $x->[$kx] - ($x_est + $self->dc);
196 2000         5557 $self->current_error($error);
197 2000         2240 my $h_new = $h;
198 2000         2356 my $tmp_coef = 1;
199 2000 50       5717 if($self->stddev_mode == 1){
200 2000         7482 $tmp_coef = $self->mu * $error / (1 + $self->stddev);
201             }else{
202 0         0 $tmp_coef = $self->mu * $error;
203             }
204 2000 50       6264 if($self->mu_mode == 1){
205 0         0 $tmp_coef = 10 * $self->mu / (1 + $self->h_length);
206             }
207              
208 2000   100     2672 for(my $k = 0;$k <= $#{$h} and $k <= $self->x_count; $k++){
  512898         2332151  
209 510898         1273502 $h_new->[$k] =
210             $h->[$k]
211             + $tmp_coef * ($x_stack->[$k] - $self->dc);
212             }
213 2000         24602 $self->h($h_new);
214             }
215             }
216              
217             ## DC component calculation and update
218             # using x_stack
219              
220             sub dc_update{
221 1000     1000 0 1492 my $self = shift;
222 1000         2299 my $x_stack = $self->x_stack;
223 1000         986 my $mean = 0;
224 1000         1383 my $num = $#$x_stack + 1;
225 1000         14812 for(0 .. $#$x_stack){
226 300000         330082 $mean += $x_stack->[$_];
227             }
228 1000         2256 $mean = $mean / $num;
229 1000         3619 $self->dc($mean);
230             }
231              
232             ## update by IIR
233             sub dc_update_iir{
234 0     0 0 0 my $self = shift;
235 0         0 my $x = shift;
236 0         0 my $new_dc = ($x - $self->dc * $self->iir_a)/(1 - $self->iir_a);
237 0         0 $self->dc($new_dc);
238             }
239              
240              
241             ## standard deviation calculation and update
242             sub stddev_update{
243 2000     2000 0 3356 my $self = shift;
244 2000         4267 my $x_stack = $self->x_stack;
245 2000         2590 my $variance = 0;
246 2000         2966 my $num = $#$x_stack + 1;
247 2000         5027 for(0 .. $#$x_stack){
248 600000         1169290 $variance += ($x_stack->[$_] - $self->dc)**2;
249             }
250 2000         3694 $variance = $variance / $num;
251 2000         7882 $self->stddev(sqrt($variance));
252             }
253              
254             ## update by IIR
255             sub stddev_update_iir{
256 0     0 0   my $self = shift;
257 0           my $x = shift;
258 0           my $tmp = sqrt(($x - $self->dc)**2);
259 0           my $new_stddev = ($tmp - $self->stddev * $self->iir_a)/(1 - $self->iir_a);
260 0           $self->stddev($new_stddev);
261             }
262              
263             ## calculation of mean value of filter
264             sub filter_dc{
265 0     0 1   my $self = shift;
266 0           my $h = $self->h;
267 0           my $mean = 0;
268 0           my $num = $#$h + 1;
269 0           for(0 .. $#$h){
270 0           $mean += $h->[$_];
271             }
272 0           return($mean / $num);
273             }
274              
275             ## calculation of stddev of filter
276             sub filter_stddev{
277 0     0 1   my $self = shift;
278 0           my $h = $self->h;
279 0           my $variance = 0;
280 0           my $num = $#$h + 1;
281 0           for(0 .. $#$h){
282 0           $variance += ($h->[$_])**2;
283             }
284 0           return(sqrt($variance / $num));
285             }
286              
287              
288              
289             1;
290             __END__