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__ |