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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
# Define some memory management.
define A Memd[a+((($2)-1)*(m+1))+($1)-1]
define B Memd[b+($1)-1]
#---------------------------------------------------------------------------
.help savgol Jun93 source
.ih
NAME
savgol -- Create a kernel for Savitzky-Golay smoothing.
.ih
USAGE
call savgol (c, np, nl, nr, ld, m)
.ih
ARGUMENTS
.ls c (O: double[np])
The smoothing kernel in "wrap-around" order. See discussion for
details.
.le
.ls np (I: int)
The number of points allocated in the array represented by the "c"
argument.
.le
.ls nl (I: int)
Size of the kernel "to the left" of the central point. See discussion
for more details.
.le
.ls nr (I: int)
Size of the kernel "to the right" of the central point. See discussion
for more details.
.le
.ls ld (I: int)
Order of the derivative desired. Should be 0 for smoothing, higher
for smoothed versions of the specified derivative.
.le
.ls m (I: int)
Order of the smoothing polynomial. Should be 0 or 1 for standard
"boxcar" or "moving window" averaging.
.le
.ih
DISCUSSION
For an introduction to Savitzky-Golay filtering, see:
.nf
Press, Teukolsky, Vetterling, & Falnnery, "Numeric Recipies:
The Art of Scientifitc Computing, Second Edition", Cambridge,
1992.
.fi
This routine returns the set of Savitzky-Golay smoothing coefficients
given the size, order of smoothing polynomial, and derivative to
return. The coefficients are returned in "wrap-around" order. Thus,
if the smoothing coefficients are C[-nl]...C[0]...C[nr], they are
returned in the array, c[i], as follows:
.nf
c[1], c[2], c[3], ..., c[nl+1],c[nl+2],...,c[np-1],c[np]
C[0], C[-1], C[-2],..., C[-nl], C[nr], ...,C[2], C[1]
.fi
A code fragment to transform the array c[i] to the orginal order,
k[i], is:
.nf
do i = 1, nl+1
k[i] = c[nl+2-i]
do i = 1, nr
k[nl+1+i] = c[np+1-i]
.fi
Array k[i], is now suitable for routines such as the IRAF VOPS routine
acnvrd.
.endhelp
#---------------------------------------------------------------------------
procedure savgol (c, np, nl, nr, ld, m)
double c[np] # O: The kernel.
int np # I: Size of the smoothing kernel.
int nl # I: Points to the left of center.
int nr # I: Points to the right of center.
int ld # I: Order of derivative to return.
int m # I: Order of the smoothing polynomial.
int imj, ipj, j, k, kk, mm, ix
double d, fac, sum
pointer indx, a, b, sp
int shifti()
begin
call smark (sp)
# Check input parameters.
if (np < nl+nr+1 || nl < 0 || nr < 0 || ld > m || nl+nr < m)
call error (1, "savgol: invalid inputs")
# Allocate memory.
call salloc (indx, m+1, TY_INT)
call salloc (a, (m+1)**2, TY_DOUBLE)
call salloc (b, m+1, TY_DOUBLE)
# Do it.
ipj = shifti (m, 1)
do ipj = 0, shifti (m, 1) {
if (ipj != 0)
sum = 0.d0
else
sum = 1.d0
do k = 1, nr
sum = sum + k**ipj
do k = 1, nl
sum = sum + (-k)**ipj
mm = min (ipj, 2*m-ipj)
do imj = -mm, mm, 2
A(1+(ipj-imj)/2,1+(ipj+imj)/2) = sum
}
call ludcmd (Memd[a], m+1, m+1, Memi[indx], d, ix)
if (ix != OK)
call error (1, "savgol: singular matrix")
do j = 1, m+1
B(j) = 0.d0
B(ld+1) = 1.d0;
call lubksd (Memd[a], m+1, m+1, Memi[indx], Memd[b])
do kk = 1, np
c[kk] = 0.d0
do k = -nl, nr {
sum = B(1)
fac = 1.d0
do mm = 1, m {
fac = fac * k
sum = sum + B(mm+1) * fac
}
kk = mod (np - k, np) + 1
c[kk] = sum
}
# That's all folks.
call sfree (sp)
end
#---------------------------------------------------------------------------
# End of savgol
#---------------------------------------------------------------------------
|
