File Coverage

blib/lib/Math/DWT.pm
Criterion Covered Total %
statement 8 168 4.7
branch 0 44 0.0
condition 0 6 0.0
subroutine 3 12 25.0
pod 9 9 100.0
total 20 239 8.3


line stmt bran cond sub pod time code
1             package Math::DWT;
2              
3 1     1   12690 use 5.006;
  1         2  
4 1     1   3 use strict;
  1         1  
  1         14  
5 1     1   3 use warnings;
  1         3  
  1         1041  
6              
7              
8             =head1 NAME
9              
10             Math::DWT - Pure Perl 1-D Discrete Wavelet Transform.
11              
12             =head1 VERSION
13              
14             Version 0.02
15              
16             =cut
17              
18             our $VERSION = '0.02';
19              
20              
21             =head1 SYNOPSIS
22              
23             This module computes the forward and inverse Discrete Wavelet Transform using arbitrary wavelets defined in modules (many are included in the distribution).
24              
25             use Math::DWT;
26              
27             my $dwt = Math::DWT->new('Daubechies',1); # Haar
28              
29             my @x = qw/8 6 7 5 3 0 9 0/; # zero-pad input to make it 2^n long
30              
31             # perform a single transform with default loaded filters
32             my $coeffs = $dwt->dwt(\@x);
33              
34             # $coeffs points to LEVEL+1 array refs
35              
36             # Check perfect reconstruction:
37             my @X = $dwt->idwt($coeffs);
38              
39             my $maxerr=0;
40             for(my $i=0;$i
41             my $err=abs($x[$i] - $X[$i]);
42             $maxerr = $err if ($err > $maxerr);
43             }
44              
45             print $maxerr . "\n"; # 5.27844434827784e-12 (close enough to 0)
46              
47              
48              
49             =head1 SUBROUTINES/METHODS
50              
51              
52             =head2 new(WAVELET,VAR)
53              
54             =head2 new(WAVELET)
55              
56             =head2 new()
57              
58             Create a new DWT object with the wavelet WAVELET and variable VAR. The wavelet is loaded from the module Math::DWT::Wavelet::WAVELET, so there will be an error if this module is not installed. If VAR is omitted, then VAR is set to whatever the default is as defined in the module (usually the lowest option). If WAVELET is left blank, then an empty object is returned with no filters defined. See the set_filters method to set custom filters (TODO).
59              
60             =cut
61              
62             sub new {
63 0     0 1   my $class=shift;
64 0           my $self={};
65 0           my $wavelet=shift;
66 0           my $var=shift;
67              
68 0           my $wname; my $w;
69 0           my ($lo_d,$hi_d,$lo_r,$hi_r);
70 0 0         if (defined($wavelet)) {
71              
72 0           $wname="Math::DWT::Wavelet::$wavelet";
73 0           require "Math/DWT/Wavelet/$wavelet.pm";
74 0           $w=$wname->new($var);
75              
76 0           ($lo_d,$hi_d,$lo_r,$hi_r)=@{$w->filters()};
  0            
77             }
78              
79 0           $self={ lo_d=>$lo_d,
80             hi_d=>$hi_d,
81             lo_r=>$lo_r,
82             hi_r=>$hi_r };
83              
84 0           bless $self,$class;
85 0           return $self;
86             };
87              
88            
89              
90             =head2 dwt(SIGNAL,LEVEL,LO-PASS,HI-PASS)
91              
92             =head2 dwt(SIGNAL,LEVEL)
93              
94             =head2 dwt(SIGNAL)
95              
96              
97             This performs the forward Discrete Wavelet Transform on SIGNAL using LO-PASS and HI-PASS as filters. The process is repeated LEVEL times. If LO-PASS and HI-PASS are omitted, then it uses defaults set in loaded wavelet. If LEVEL is omitted, then it sets LEVEL=1. SIGNAL is an array ref. The length of SIGNAL must be a power of 2, e.g. 256, 512, 1024, etc.
98              
99             dwt returns an arrayref (or array, if called in array context) of arrayrefs, one for each set of detail coefficients, and an extra one for the last set of scaling coefficients. So for a 256-sample long input signal, processed to 4 levels, the return value of dwt would be an array/arrayref with 5 arrayrefs with sizes 128, 64, 32, 16, and 16, in that order.
100              
101             @x=@x[0..255];
102             my $coeffs = $dwt->dwt(\@x,4);
103              
104             foreach my $i (@$coeffs) { print scalar(@$i) . " "; }
105             # prints
106             # 128 64 32 16 16
107              
108              
109              
110             =cut
111              
112             sub dwt {
113 0     0 1   my $self=shift;
114 0           my ($x,$level,$lo,$hi)=@_;
115              
116 0 0 0       if (!defined($level)|| $level == 0) {
117 0           $level = 1;
118             };
119              
120 0 0         if (!defined($lo)) {
121 0           $lo = $self->{lo_d};
122             };
123 0 0         if (!defined($hi)) {
124 0           $hi = $self->{hi_d};
125             };
126              
127 0           my @hi_tmp;
128 0           my @tmp_lo=@{$x};
  0            
129 0           my $lo_tmp=\@tmp_lo;
130              
131 0           for(my $j=1;$j<=$level;$j++) {
132 0           ($lo_tmp,$hi_tmp[$j-1])=afb($lo_tmp,$lo,$hi);
133             };
134 0           $hi_tmp[$level]=$lo_tmp;
135 0 0         if (wantarray()) {
136 0           return @hi_tmp;
137             };
138 0           return \@hi_tmp;
139              
140             }
141              
142             =head2 idwt(COEFFICIENTS,LEVEL,LO-PASS,HI-PASS)
143              
144             =head2 idwt(COEFFICIENTS,LEVEL)
145              
146             =head2 idwt(COEFFICIENTS)
147              
148             Inverse Discrete Wavelet Transform. Same defaults setup as for dwt. COEFFICIENTS must be an array ref of the structure returned by dwt. Returns an array or arrayref (depending on the value of wantarray()) with the original signal, or at least the signal inversely processed LEVEL times.
149              
150             =cut
151              
152             sub idwt {
153 0     0 1   my $self=shift;
154 0           my ($x,$level,$lo,$hi)=@_;
155              
156 0 0 0       if (!defined($level)|| $level == 0) {
157 0           $level = 1;
158             };
159              
160 0 0         if (!defined($lo)) {
161 0           $lo = $self->{lo_r};
162             };
163 0 0         if (!defined($hi)) {
164 0           $hi = $self->{hi_r};
165             };
166              
167 0           my @y;
168              
169 0           @y=@{$x->[$level]};
  0            
170             #print STDERR "HERE: " . scalar(@y) . "\n";
171              
172 0           for (my $j=$level;$j>=1;$j--) {
173 0           @y=sfb(\@y,$x->[$j-1],$lo,$hi);
174             };
175              
176 0           my @out;
177 0           for (my $i=0;$i
178 0           $out[$i]=$y[$i];
179             };
180 0 0         if (wantarray()) {
181 0           return (@out);
182             };
183 0           return \@out;
184            
185             }
186              
187             =head2 cshift(ARRAY,SHIFT)
188              
189             Utility method that shifts the elements of ARRAY by SHIFT amount in a circular fashion (i.e. elements positioned past the end of the array are placed back at the beginning). Accepts negative values for SHIFT.
190              
191             =cut
192              
193              
194             sub cshift {
195 0     0 1   my ($xc,$m)=@_;
196 0           my @ac=@{$xc};
  0            
197 0           my @bc;
198             my @tmpc;
199 0 0         if ($m < 0) {
    0          
200 0           @tmpc=splice(@ac,0,-$m);
201 0           @bc=(@ac,@tmpc);
202             } elsif ($m > 0) {
203 0           @tmpc=splice(@ac,-$m);
204 0           @bc=(@tmpc,@ac);
205             } else {
206 0           @bc=@ac;
207             }
208              
209 0 0         if (wantarray()) {
210 0           return @bc;
211             };
212 0           return \@bc;
213             };
214              
215             =head2 afb(SIGNAL,LO-PASS,HI-PASS)
216              
217             This is a utility method you shouldn't need to use. Stands for Analysis Filter Bank. This method convolves SIGNAL with LO-PASS and separately with HI-PASS, then downsamples each by 2, and returns arrayrefs of the results of the proper length (half the input length).
218              
219             =cut
220              
221             sub afb {
222             # analysis filter banks
223 0     0 1   my ($x,$afL,$afH)=@_;
224 0           my $N=scalar(@{$x});
  0            
225 0           my $L=scalar(@{$afL})/2;
  0            
226              
227             # do the mystery circle shift
228 0           my @X=cshift($x,-$L);
229             #foreach(@X) { $_=$_/sqrt(2); }
230              
231             # temporary lo/hi arrays
232 0           my @c;
233             my @d;
234              
235             # convolutions
236 0           for (my $i=0; $i < $N; $i++) {
237 0           $c[$i]=0;
238 0           $d[$i]=0;
239 0           for (my $j = 0; $j < $L*2; $j++) {
240 0           my $a=0;
241 0 0         if (defined($X[$i-$j])) { $a=$X[$i-$j]; }
  0            
242 0           $c[$i]+=$a * $afL->[$j];
243 0           $d[$i]+=$a * $afH->[$j];
244             };
245             };
246              
247 0           my @lo; my @hi;
248 0           foreach my $k (0 .. scalar(@c)/2) {
249 0           push @lo, $c[2*$k];
250 0           push @hi, $d[2*$k];
251             };
252 0           for (my $i=0;$i<$L;$i++) {
253 0 0         if (!defined($lo[$i+($N/2)])) { $lo[$i+($N/2)]=0;}
  0            
254 0 0         if (!defined($hi[$i+($N/2)])) { $hi[$i+($N/2)]=0;}
  0            
255 0           $lo[$i]=$lo[$i+($N/2)] + $lo[$i];
256 0           $hi[$i]=$hi[$i+($N/2)] + $hi[$i];
257             }
258              
259 0           my @lo_fin=@lo[ 0 .. ($N/2)-1 ];
260 0           my @hi_fin=@hi[ 0 .. ($N/2)-1 ];
261              
262 0           return (\@lo_fin, \@hi_fin);
263             };
264              
265             =head2 sfb(LO-PART,HI-PART,LO-PASS,HI-PASS)
266              
267             This is an internal utility method for reconstructing a signal from a set of wavelet coefficients and scaling coefficients and a pair of lo-pass and hi-pass reconstruction filters. It stands for Synthesis Filter Bank and it returns a single array or arrayref that is the result of upsampling the input by 2, then convolving each set of coefficients with its respective filter, and, finally, adding up corresponding lo-pass/hi-pass values.
268              
269             =cut
270              
271              
272             sub sfb {
273 0     0 1   my ($lo,$hi,$sfL,$sfH)=@_;
274 0           my $N=2*scalar(@{$lo});
  0            
275 0           my $L=scalar(@{$sfL});
  0            
276 0           $lo=upfirdn($lo,$sfL,2,1);
277 0           $hi=upfirdn($hi,$sfH,2,1);
278 0           my @y;
279 0           my $i=0;
280              
281 0           for ($i=0;$i
  0            
282 0           $y[$i]=$lo->[$i] + $hi->[$i];
283             };
284              
285 0           for ($i=0;$i<$L-$N/2;$i++) {
286 0           my $a=0;
287 0 0         if (defined($y[$N+$i])) { $a=$y[$N+$i]; }
  0            
288 0           my $b=0;
289 0 0         if (defined($y[$i])) { $b=$y[$i]; }
  0            
290 0           $y[$i]=$b + $a;
291             };
292 0           @y=@y[0..$N-1];
293 0           @y=cshift(\@y,1-($L/2));
294             #foreach(@y) { $_=$_/sqrt(2); }
295              
296 0 0         if (wantarray()) { return @y; };
  0            
297 0           return \@y;
298             };
299              
300             =head2 upfirdn(A,B,UP,DN)
301              
302             This is an internal utility method modeled after MATLAB's upfirdn command. First, it upsamples signal A by the value of UP (use 1 to disable), then convolves A and B, then downsamples the result by a factor of DN (again, use 1 to disable). The result is an array or arrayref with length defined by: (((length(A)-1)*UP+length(B))/DN)
303              
304             =cut
305              
306              
307              
308             sub upfirdn {
309 0     0 1   my ($a,$b,$up,$dn)=@_;
310 0           my $N=scalar(@{$a});
  0            
311 0           my $L=scalar(@{$b});
  0            
312              
313 0           my @A;
314 0 0         if ($up > 1) {
315 0           for(my $k=0; $k
  0            
316 0           $A[$up*$k]=$a->[$k];
317 0           $A[($up*$k)+1]=0;
318             };
319             } else {
320 0           for(my $k=0;$k
  0            
321 0           $A[$k]=$a->[$k];
322             };
323             };
324              
325 0           my @c;
326 0           for (my $i=0; $i < scalar(@A); $i++) {
327 0           $c[$i]=0;
328 0           for (my $j = 0; $j < $L; $j++) {
329 0           my $tmpa=0;
330 0 0         if (defined($A[$i-$j])) { $tmpa=$A[$i-$j]; }
  0            
331 0           $c[$i]+=$tmpa * $b->[$j];
332             };
333             };
334 0           my @d;
335 0           foreach my $k (0 .. (scalar(@c)/$dn)-1) {
336 0           push @d, $c[$dn*$k];
337             };
338 0 0         if (wantarray()) {
339 0           return @d;
340             }
341 0           return \@d;
342             };
343              
344              
345             =head2 set_filters(LO_D,HI_D,LO_R,HI_R)
346              
347             Set the filters manually. Each of LO_D, HI_D, LO_R, and HI_R is an arrayref to the set of coefficients for the respective filter set. Returns undef.
348              
349             =cut
350              
351             sub set_filters {
352 0     0 1   my $self=shift;
353 0           my ($lo_d,$hi_d,$lo_r,$hi_r)=@_;
354              
355 0           $self->{lo_d}=$lo_d;
356 0           $self->{hi_d}=$hi_d;
357 0           $self->{lo_r}=$lo_r;
358 0           $self->{hi_r}=$hi_r;
359              
360 0           return undef;
361             };
362              
363             =head2 get_filters()
364              
365             Get the set of filters. Returns an array or an arrayref containing LO_D, HI_D, LO_R, HI_R in that order.
366              
367             =cut
368              
369             sub get_filters {
370 0     0 1   my $self=shift;
371 0           my ($lo_d,$hi_d,$lo_r,$hi_r)=($self->{lo_d},$self->{hi_d},$self->{lo_r},$self->{hi_r});
372              
373 0 0         if (wantarray()) {
374 0           return ($lo_d,$hi_d,$lo_r,$hi_r);
375             };
376 0           return [$lo_d, $hi_d, $lo_r, $hi_r ];
377             };
378              
379              
380              
381              
382              
383              
384             =head1 SEE ALSO
385              
386             Math::DWT::UDWT(3pm), Math::DWT::Wavelet::Haar(3pm), Math::DWT::Wavelet::Coiflet(3pm), Math::DWT::Wavelet::Symlet(3pm), Math::DWT::Wavelet::Biorthogonal(3pm), Math::DWT::Wavelet::ReverseBiorthogonal(3pm), Math::DWT::Wavelet::Daubechies(3pm), Math::DWT::Wavelet::DiscreteMeyer(3pm), perl(1)
387              
388              
389              
390             =head1 AUTHOR
391              
392              
393             Mike Kroh, C<< >>
394              
395             =head1 BUGS
396              
397              
398             Please report any bugs or feature requests to C, or through
399             the web interface at L. I will be notified, and then you'll
400             automatically be notified of progress on your bug as I make changes.
401              
402              
403              
404             =head1 ACKNOWLEDGEMENTS
405              
406              
407             Special thanks to Ivan Selesnick for his software (on which these modules are based) available here L.
408              
409             Some wavelet filter coefficients scraped from this site: L.
410              
411             =head1 LICENSE AND COPYRIGHT
412              
413              
414             Copyright 2016 Mike Kroh.
415              
416             This program is free software; you can redistribute it and/or modify it
417             under the terms of the the Artistic License (2.0). You may obtain a
418             copy of the full license at:
419              
420             L
421              
422             Any use, modification, and distribution of the Standard or Modified
423             Versions is governed by this Artistic License. By using, modifying or
424             distributing the Package, you accept this license. Do not use, modify,
425             or distribute the Package, if you do not accept this license.
426              
427             If your Modified Version has been derived from a Modified Version made
428             by someone other than you, you are nevertheless required to ensure that
429             your Modified Version complies with the requirements of this license.
430              
431             This license does not grant you the right to use any trademark, service
432             mark, tradename, or logo of the Copyright Holder.
433              
434             This license includes the non-exclusive, worldwide, free-of-charge
435             patent license to make, have made, use, offer to sell, sell, import and
436             otherwise transfer the Package with respect to any patent claims
437             licensable by the Copyright Holder that are necessarily infringed by the
438             Package. If you institute patent litigation (including a cross-claim or
439             counterclaim) against any party alleging that the Package constitutes
440             direct or contributory patent infringement, then this Artistic License
441             to you shall terminate on the date that such litigation is filed.
442              
443             Disclaimer of Warranty: THE PACKAGE IS PROVIDED BY THE COPYRIGHT HOLDER
444             AND CONTRIBUTORS "AS IS' AND WITHOUT ANY EXPRESS OR IMPLIED WARRANTIES.
445             THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
446             PURPOSE, OR NON-INFRINGEMENT ARE DISCLAIMED TO THE EXTENT PERMITTED BY
447             YOUR LOCAL LAW. UNLESS REQUIRED BY LAW, NO COPYRIGHT HOLDER OR
448             CONTRIBUTOR WILL BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, OR
449             CONSEQUENTIAL DAMAGES ARISING IN ANY WAY OUT OF THE USE OF THE PACKAGE,
450             EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
451              
452              
453             =cut
454              
455             1; # End of Math::DWT