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
|
#---------------------------------------------------------------------------
.help sg_convolve Jun93 source
.ih
NAME
sg_convolve -- Convolve an array using Savitzky-Golay filter.
.ih
USAGE
call sg_convolve (size, order, in, out, n)
.ih
ARGUMENTS
.ls size (I: int)
The full size of the smoothing kernel. If less than or equal to 1,
then no convolving takes place.
.le
.ls order (I: int)
The order of the smoothing polynomial. For normal "boxcar" smoothing,
this should be 0 or 1. Greater values preserve higher order terms in the
original data. Larger sizes are needed for this to be effective.
.le
.ls in (I: double[n])
The data array to be convolved.
.le
.ls out (O: double[n])
The convolved array. May be the same as the input array.
.le
.ih
DESCRIPTION
The routine, savgol, is used to calculate a Savitsky-Golay convolving
kernel. This kernel is then applied, using standard routines, to the
input data. See the routine savgol for more information.
.ih
SEE ALSO
savgol
.endhelp
#---------------------------------------------------------------------------
procedure sg_convolve (size, order, in, out, n)
int size # I: The size of the filter.
int order # I: The order to preserve while filtering.
double in[n] # I: Data to be convolved.
double out[n] # O: The convolved data.
int n # I: Length of the arrays.
# Kernel parameters.
int half # Half size of kernel.
int isize # Odd size of kernel.
pointer k, kx # The kernel in real/double-wrap versions.
# Misc.
pointer adx # Generic double array.
int i # Generic.
pointer sp # Stack pointer.
begin
call smark (sp)
# Fix the kernel size to be odd.
half = size / 2
isize = half * 2 + 1
# Make sure there is something to convolve. If not, just copy and
# run.
if (isize <= 1)
call amovd (in, out, n)
else {
call salloc (k, isize, TY_REAL)
call salloc (kx, isize, TY_DOUBLE)
call salloc (adx, n+isize, TY_DOUBLE)
# Compute the kernel.
call savgol (Memd[kx], isize, half, half, 0, order)
do i = 0, half
Memr[k+i] = Memd[kx+half-i]
do i = 0, half-1
Memr[k+half+i+1] = Memd[kx+isize -i-1]
# Put the data in the extended array and pad the ends as
# constants.
call amovd (in, Memd[adx+half], n)
do i = 1, half {
Memd[adx+half-i] = in[1]
Memd[adx+half+n+i-1] = in[n]
}
# Filter it.
call acnvrd (Memd[adx], out, n, Memr[k], isize)
}
# That's all folks.
call sfree (sp)
end
#---------------------------------------------------------------------------
# End of sg_convolve
#---------------------------------------------------------------------------
|