1 |
chrisfen |
1287 |
/*************************************************************************** |
2 |
|
|
************************************************************************** |
3 |
|
|
|
4 |
|
|
S2kit 1.0 |
5 |
|
|
|
6 |
|
|
A lite version of Spherical Harmonic Transform Kit |
7 |
|
|
|
8 |
|
|
Peter Kostelec, Dan Rockmore |
9 |
|
|
{geelong,rockmore}@cs.dartmouth.edu |
10 |
|
|
|
11 |
|
|
Contact: Peter Kostelec |
12 |
|
|
geelong@cs.dartmouth.edu |
13 |
|
|
|
14 |
|
|
Copyright 2004 Peter Kostelec, Dan Rockmore |
15 |
|
|
|
16 |
|
|
This file is part of S2kit. |
17 |
|
|
|
18 |
|
|
S2kit is free software; you can redistribute it and/or modify |
19 |
|
|
it under the terms of the GNU General Public License as published by |
20 |
|
|
the Free Software Foundation; either version 2 of the License, or |
21 |
|
|
(at your option) any later version. |
22 |
|
|
|
23 |
|
|
S2kit is distributed in the hope that it will be useful, |
24 |
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of |
25 |
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
26 |
|
|
GNU General Public License for more details. |
27 |
|
|
|
28 |
|
|
You should have received a copy of the GNU General Public License |
29 |
|
|
along with S2kit; if not, write to the Free Software |
30 |
|
|
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
31 |
|
|
|
32 |
|
|
See the accompanying LICENSE file for details. |
33 |
|
|
|
34 |
|
|
************************************************************************ |
35 |
|
|
************************************************************************/ |
36 |
|
|
|
37 |
|
|
|
38 |
|
|
/*************************************************************************/ |
39 |
|
|
|
40 |
|
|
/* Source code to synthesize functions using a naive method |
41 |
|
|
based on recurrence. This is slow but does not require any |
42 |
|
|
precomputed functions, and is also stable. |
43 |
|
|
*/ |
44 |
|
|
|
45 |
|
|
|
46 |
|
|
#include <math.h> |
47 |
|
|
#include <string.h> |
48 |
|
|
|
49 |
|
|
|
50 |
|
|
/************************************************************************/ |
51 |
|
|
/* |
52 |
|
|
|
53 |
|
|
Naive_AnalysisX: computing the discrete Legendre transform of |
54 |
|
|
a function via summing naively. I.e. This is |
55 |
|
|
the FORWARD discrete Legendre transform. |
56 |
|
|
|
57 |
|
|
bw - bandwidth |
58 |
|
|
m - order |
59 |
|
|
|
60 |
|
|
data - a pointer to double array of size (2*bw) containing |
61 |
|
|
the sample points |
62 |
|
|
|
63 |
|
|
result - a pointer to double array of size (bw-m) which, at the |
64 |
|
|
conclusion of the routine, will contains the coefficients |
65 |
|
|
|
66 |
|
|
plmtable - a pointer to a double array of size (2*bw*(bw-m)); |
67 |
|
|
contains the PRECOMPUTED plms, i.e. associated Legendre |
68 |
|
|
functions. E.g. Should be generated by a call to |
69 |
|
|
|
70 |
|
|
PmlTableGen() |
71 |
|
|
|
72 |
|
|
(see pmls.c) |
73 |
|
|
|
74 |
|
|
NOTE that these Legendres are normalized with norm |
75 |
|
|
equal to 1 !!! |
76 |
|
|
|
77 |
|
|
workspace - array of size 2 * bw; |
78 |
|
|
|
79 |
|
|
|
80 |
|
|
*/ |
81 |
|
|
|
82 |
|
|
|
83 |
|
|
void Naive_AnalysisX(double *data, |
84 |
|
|
int bw, |
85 |
|
|
int m, |
86 |
|
|
double *weights, |
87 |
|
|
double *result, |
88 |
|
|
double *plmtable, |
89 |
|
|
double *workspace) |
90 |
|
|
{ |
91 |
|
|
int i, j; |
92 |
|
|
double result0, result1, result2, result3; |
93 |
|
|
register double *wdata; |
94 |
|
|
|
95 |
|
|
wdata = workspace; |
96 |
|
|
|
97 |
|
|
/* make sure result is zeroed out */ |
98 |
|
|
memset( result, 0, sizeof(double) * (bw - m) ); |
99 |
|
|
|
100 |
|
|
/* apply quadrature weights */ |
101 |
|
|
/* |
102 |
|
|
I only have to differentiate between even and odd |
103 |
|
|
weights when doing something like seminaive, something |
104 |
|
|
which involves the dct. In this naive case, the parity of |
105 |
|
|
the order of the transform doesn't matter because I'm not |
106 |
|
|
dividing by sin(x) when precomputing the Legendres (because |
107 |
|
|
I'm not taking their dct). The plain ol' weights are just |
108 |
|
|
fine. |
109 |
|
|
*/ |
110 |
|
|
|
111 |
|
|
for(i = 0; i < 2 * bw; i++) |
112 |
|
|
wdata[i] = data[i] * weights[i]; |
113 |
|
|
|
114 |
|
|
/* unrolling seems to work */ |
115 |
|
|
if ( 1 ) |
116 |
|
|
{ |
117 |
|
|
for (i = 0; i < bw - m; i++) |
118 |
|
|
{ |
119 |
|
|
result0 = 0.0; result1 = 0.0; |
120 |
|
|
result2 = 0.0; result3 = 0.0; |
121 |
|
|
|
122 |
|
|
for(j = 0; j < (2 * bw) % 4; ++j) |
123 |
|
|
result0 += wdata[j] * plmtable[j]; |
124 |
|
|
|
125 |
|
|
for( ; j < (2 * bw); j += 4) |
126 |
|
|
{ |
127 |
|
|
result0 += wdata[j] * plmtable[j]; |
128 |
|
|
result1 += wdata[j + 1] * plmtable[j + 1]; |
129 |
|
|
result2 += wdata[j + 2] * plmtable[j + 2]; |
130 |
|
|
result3 += wdata[j + 3] * plmtable[j + 3]; |
131 |
|
|
} |
132 |
|
|
result[i] = result0 + result1 + result2 + result3; |
133 |
|
|
|
134 |
|
|
plmtable += (2 * bw); |
135 |
|
|
} |
136 |
|
|
} |
137 |
|
|
else |
138 |
|
|
{ |
139 |
|
|
for (i = 0; i < bw - m; i++) |
140 |
|
|
{ |
141 |
|
|
result0 = 0.0 ; |
142 |
|
|
for(j = 0; j < (2 * bw) ; j++) |
143 |
|
|
result0 += wdata[j] * plmtable[j]; |
144 |
|
|
result[i] = result0 ; |
145 |
|
|
|
146 |
|
|
plmtable += (2 * bw); |
147 |
|
|
} |
148 |
|
|
} |
149 |
|
|
} |
150 |
|
|
|
151 |
|
|
|
152 |
|
|
/************************************************************************/ |
153 |
|
|
/* This is the procedure that synthesizes a function from a list |
154 |
|
|
of coefficients of a Legendre series. I.e. this is the INVERSE |
155 |
|
|
discrete Legendre transform. |
156 |
|
|
|
157 |
|
|
Function is synthesized at the (2*bw) Chebyshev nodes. Associated |
158 |
|
|
Legendre functions are assumed to be precomputed. |
159 |
|
|
|
160 |
|
|
bw - bandwidth |
161 |
|
|
|
162 |
|
|
m - order |
163 |
|
|
|
164 |
|
|
coeffs - a pointer to double array of size (bw-m). First coefficient is |
165 |
|
|
coefficient for Pmm |
166 |
|
|
result - a pointer to double array of size (2*bw); at the conclusion |
167 |
|
|
of the routine, this array will contain the |
168 |
|
|
synthesized function |
169 |
|
|
|
170 |
|
|
plmtable - a pointer to a double array of size (2*bw*(bw-m)); |
171 |
|
|
contains the PRECOMPUTED plms, i.e. associated Legendre |
172 |
|
|
functions. E.g. Should be generated by a call to |
173 |
|
|
|
174 |
|
|
PmlTableGen(), |
175 |
|
|
|
176 |
|
|
(see pmls.c) |
177 |
|
|
|
178 |
|
|
|
179 |
|
|
NOTE that these Legendres are normalized with norm |
180 |
|
|
equal to 1 !!! |
181 |
|
|
*/ |
182 |
|
|
|
183 |
|
|
void Naive_SynthesizeX(double *coeffs, |
184 |
|
|
int bw, |
185 |
|
|
int m, |
186 |
|
|
double *result, |
187 |
|
|
double *plmtable) |
188 |
|
|
{ |
189 |
|
|
int i, j; |
190 |
|
|
double tmpcoef; |
191 |
|
|
|
192 |
|
|
/* make sure result is zeroed out */ |
193 |
|
|
memset( result, 0, sizeof(double) * 2 * bw ); |
194 |
|
|
|
195 |
|
|
for ( i = 0 ; i < bw - m ; i ++ ) |
196 |
|
|
{ |
197 |
|
|
tmpcoef = coeffs[i] ; |
198 |
|
|
if ( tmpcoef != 0.0 ) |
199 |
|
|
for (j=0; j<(2*bw); j++) |
200 |
|
|
result[j] += (tmpcoef * plmtable[j]); |
201 |
|
|
plmtable += (2 * bw ) ; |
202 |
|
|
} |
203 |
|
|
} |