1 |
/*************************************************************************** |
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 |
} |