|
1 #====================================================================== |
|
2 # . I N T E R P . A D C P |
|
3 # doc: Fri Apr 16 16:07:48 2010 |
|
4 # dlm: Tue Aug 9 23:07:35 2011 |
|
5 # (c) 2010 A.M. Thurnherr |
|
6 # uE-Info: 81 0 NIL 0 0 72 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # interpolation scheme mimicking RDI ADCP response as described in |
|
10 # Broadband primer, p.17 |
|
11 |
|
12 # HISTORY: |
|
13 # Apr 16, 2010: - created |
|
14 # Aug 9, 2011: - added -u |
|
15 |
|
16 # NOTES: |
|
17 # - interface is described in [.interp.linear] |
|
18 |
|
19 $IS_opts = 'b:u'; |
|
20 $IS_optsUsage = '[pass -u)nfiltered] -b)in/pulse <length>'; |
|
21 |
|
22 sub IS_usage() |
|
23 { |
|
24 &antsUsageError() |
|
25 unless defined($opt_b); |
|
26 } |
|
27 |
|
28 sub IS_init($$$$) {} |
|
29 |
|
30 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
31 # |
|
32 # &IS_interpolate(br,idr,xf,xv,xi,f) interpolate field f |
|
33 # br data buffer reference |
|
34 # idr init-data reference |
|
35 # xf x field |
|
36 # xv x value |
|
37 # xi index of last record with x-value <= x |
|
38 # f field number to interpolate |
|
39 # <ret val> interpolated value |
|
40 # |
|
41 # NB: |
|
42 # - handle f == xf |
|
43 # - return NaN if any of the y values required is NaN |
|
44 # |
|
45 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
46 |
|
47 sub IS_interpolate($$$$$$) |
|
48 { |
|
49 my($bR,$idR,$xf,$xv,$xi,$f) = @_; |
|
50 |
|
51 return $xv if ($xf == $f); # return target x |
|
52 |
|
53 my($tow) = $xi; # top of triangular sampling window |
|
54 while ($tow>=0 && (!numberp($bR->[$tow][$xf]) || $bR->[$tow][$xf] > $xv-$opt_b)) { |
|
55 $tow--; |
|
56 } |
|
57 if ($tow < 0) { # incomplete window |
|
58 return nan unless ($opt_u); |
|
59 $tow = 0; |
|
60 } else { |
|
61 $tow++; |
|
62 } |
|
63 |
|
64 my($bow) = $xi+1; # bottom of triangular sampling window |
|
65 while ($bow<=$#{$bR} && (!numberp($bR->[$bow][$xf]) || $bR->[$bow][$xf] < $xv+$opt_b)) { |
|
66 $bow++; |
|
67 } |
|
68 |
|
69 if ($bow > $#{$bR}) { # incomplete window |
|
70 return nan unless ($opt_u); |
|
71 $bow = $#{$bR}; |
|
72 } else { |
|
73 $bow--; |
|
74 } |
|
75 |
|
76 my($sweight) = 0; # calculate weighted average |
|
77 my($sum) = 0; |
|
78 for (my($i)=$tow; $i<=$bow; $i++) { |
|
79 next unless (numberp($bR->[$i][$xf]) && numberp($bR->[$i][$f])); |
|
80 my($weight) = 1 - abs($bR->[$i][$xf]-$xv)/$opt_b; |
|
81 $sum += $weight * $bR->[$i][$f]; |
|
82 $sweight += $weight; |
|
83 } |
|
84 |
|
85 return ($sweight>0) ? $sum/$sweight : nan; |
|
86 } |
|
87 |
|
88 #---------------------------------------------------------------------- |
|
89 |
|
90 1; |