SAF
Loading...
Searching...
No Matches
test__sh_module.c
Go to the documentation of this file.
1/*
2 * Copyright 2020-2021 Leo McCormack
3 *
4 * Permission to use, copy, modify, and/or distribute this software for any
5 * purpose with or without fee is hereby granted, provided that the above
6 * copyright notice and this permission notice appear in all copies.
7 *
8 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
9 * REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
10 * AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
11 * INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
12 * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
13 * OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
14 * PERFORMANCE OF THIS SOFTWARE.
15 */
16
25#include "saf_test.h" /* the SAF unit test declarations */
26
27void test__getSHreal(void){
28 int i, j, k, order, nDirs, nSH;
29 float scale;
30 float* t_dirs_deg;
31 float** t_dirs_rad, **Y, **YYT;
32
33 /* Config */
34 const float acceptedTolerance = 0.00001f;
35 int nTestOrders = 10;
36 int testOrders[10] = {1,2,3,4,5,6,7,8,9,10};
37
38 /* Loop over orders */
39 for(i=0; i<nTestOrders; i++){
40 order = testOrders[i];
41 nSH = ORDER2NSH(order);
42
43 /* Pull an appropriate t-design */
44 t_dirs_deg = (float*)__HANDLES_Tdesign_dirs_deg[2*order];
45 nDirs = __Tdesign_nPoints_per_degree[2*order];
46 t_dirs_rad = (float**)malloc2d(nDirs, 2, sizeof(float));
47 for(j=0; j<nDirs; j++){
48 t_dirs_rad[j][0] = t_dirs_deg[j*2] * SAF_PI/180.0f;
49 t_dirs_rad[j][1] = SAF_PI/2.0f - t_dirs_deg[j*2+1] * SAF_PI/180.0f; /* elevation->inclination */
50 }
51
52 /* Compute spherical harmonic coefficients */
53 Y = (float**)malloc2d(nSH, nDirs, sizeof(float));
54 getSHreal(order, FLATTEN2D(t_dirs_rad), nDirs, FLATTEN2D(Y));
55 scale = SQRT4PI;
56 utility_svsmul(FLATTEN2D(Y), &scale, nSH*nDirs, FLATTEN2D(Y));
57
58 /* Check Y is orthogonal: */
59 YYT = (float**)malloc2d(nSH, nSH, sizeof(float));
60 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, nSH, nSH, nDirs, 1.0f,
61 FLATTEN2D(Y), nDirs,
62 FLATTEN2D(Y), nDirs, 0.0f,
63 FLATTEN2D(YYT), nSH);
64
65 /* Should be Identity: */
66 scale = 1.0f/(float)nDirs;
67 utility_svsmul(FLATTEN2D(YYT), &scale, nSH*nSH, FLATTEN2D(YYT));
68 for(j=0; j<nSH; j++){
69 for(k=0; k<nSH; k++) {
70 if(j==k)
71 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 1.0f, YYT[j][k]);
72 else
73 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 0.0f, YYT[j][k]);
74 }
75 }
76
77 /* clean-up */
78 free(t_dirs_rad);
79 free(Y);
80 free(YYT);
81 }
82}
83
84void test__getSHreal_recur(void){
85 int i, j;
86
87 /* Config */
88 /* In general, the values from this recusive alternative are well below this
89 * tolerance value. However, the error does get larger for higher-orders and
90 * when dir[1] is near 0. */
91 float acceptedTolerance = 0.005f;
92 const int order = 15;
93 const int nSH = ORDER2NSH(order);
94 float dir[2];
95
96 /* Check that the output of getSHreal_recur matches that of getSH_recur */
97 float Yr[ORDER2NSH(15)];
98 float Y[ORDER2NSH(15)];
99 for(i=0; i<1e3; i++){
100 rand_m1_1(&dir[0] , 1);
101 rand_m1_1(&dir[1] , 1);
102 dir[0] *= SAF_PI;
103 dir[1] *= SAF_PI/2.0f;
104 getSHreal_recur(order, (float*)dir, 1, (float*)Yr);
105 getSHreal(order, (float*)dir, 1, (float*)Y);
106 for(j=0; j<nSH; j++)
107 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, Yr[j], Y[j]);
108 }
109}
110
111void test__getSHreal_part(void){
112 int i, j;
113
114 /* Config */
115 /* In general, the values from this recusive alternative are well below this
116 * tolerance value. However, the error does get larger for higher-orders and
117 * when dir[1] is near 0. */
118 float acceptedTolerance = 0.005f;
119 const int order_start = 5;
120 const int order_end = 15;
121 const int nSH_end = ORDER2NSH(order_end);
122 float dir[2];
123
124 /* Check that the output of getSHreal_part matches that of getSH_recur */
125 float Yp0[ORDER2NSH(15)];
126 float Yp5[ORDER2NSH(15)];
127 float Y[ORDER2NSH(15)];
128 for(i=0; i<1e3; i++){
129 rand_m1_1(&dir[0] , 1);
130 rand_m1_1(&dir[1] , 1);
131 dir[0] *= SAF_PI;
132 dir[1] *= SAF_PI/2.0f;
133 getSHreal_part(0, order_end, (float*)dir, 1, (float*)Yp0);
134 getSHreal_part(order_start, order_end, (float*)dir, 1, (float*)Yp5);
135 getSHreal(order_end, (float*)dir, 1, (float*)Y);
136 for(j=0; j<nSH_end; j++)
137 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, Yp0[j], Y[j]);
138 for(j=(order_start+1)*(order_start+1); j<nSH_end; j++)
139 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, Yp5[j], Y[j]);
140 }
141}
142
143void test__getSHcomplex(void){
144 int i, j, k, order, nDirs, nSH;
145 float_complex scale;
146 float* t_dirs_deg;
147 float** t_dirs_rad;
148 float_complex **Y, **YYH;
149 const float_complex calpha = cmplxf(1.0f, 0.0f), cbeta = cmplxf(0.0f, 0.0f);
150
151 /* Config */
152 const float acceptedTolerance = 0.00001f;
153 int nTestOrders = 10;
154 int testOrders[10] = {1,2,3,4,5,6,7,8,9,10};
155
156 /* Loop over orders */
157 for(i=0; i<nTestOrders; i++){
158 order = testOrders[i];
159 nSH = ORDER2NSH(order);
160
161 /* Pull an appropriate t-design */
162 t_dirs_deg = (float*)__HANDLES_Tdesign_dirs_deg[2*order];
163 nDirs = __Tdesign_nPoints_per_degree[2*order];
164 t_dirs_rad = (float**)malloc2d(nDirs, 2, sizeof(float));
165 for(j=0; j<nDirs; j++){
166 t_dirs_rad[j][0] = t_dirs_deg[j*2] * SAF_PI/180.0f;
167 t_dirs_rad[j][1] = SAF_PI/2.0f - t_dirs_deg[j*2+1] * SAF_PI/180.0f; /* elevation->inclination */
168 }
169
170 /* Compute spherical harmonic coefficients */
171 Y = (float_complex**)malloc2d(nSH, nDirs, sizeof(float_complex));
172 getSHcomplex(order, FLATTEN2D(t_dirs_rad), nDirs, FLATTEN2D(Y));
173 scale = cmplxf(SQRT4PI, 0.0f);
174 utility_cvsmul(FLATTEN2D(Y), &scale, nSH*nDirs, FLATTEN2D(Y));
175
176 /* Check Y is orthogonal: */
177 YYH = (float_complex**)malloc2d(nSH, nSH, sizeof(float_complex));
178 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, nSH, nSH, nDirs, &calpha,
179 FLATTEN2D(Y), nDirs,
180 FLATTEN2D(Y), nDirs, &cbeta,
181 FLATTEN2D(YYH), nSH);
182
183 /* Should be Identity: */
184 scale = cmplxf(1.0f/(float)nDirs, 0.0f);
185 utility_cvsmul(FLATTEN2D(YYH), &scale, nSH*nSH, FLATTEN2D(YYH));
186 for(j=0; j<nSH; j++){
187 for(k=0; k<nSH; k++) {
188 if(j==k)
189 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 1.0f, crealf(YYH[j][k]));
190 else
191 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 0.0f, crealf(YYH[j][k]));
192 }
193 }
194
195 /* clean-up */
196 free(t_dirs_rad);
197 free(Y);
198 free(YYH);
199 }
200}
201
202void test__getSHrotMtxReal(void){
203 int i,j,nSH,order;
204 float Rzyx[3][3];
205 float** Mrot;
206
207 /* Config */
208 const float acceptedTolerance = 0.00001f;
209
210 /* Rotation matrix for 0,0,0 should be identity */
211 order = 22;
212 yawPitchRoll2Rzyx(0.0f, 0.0f, 0.0f, 0, Rzyx);
213 nSH = ORDER2NSH(order);
214 Mrot = (float**)malloc2d(nSH, nSH, sizeof(float));
215 getSHrotMtxReal(Rzyx, FLATTEN2D(Mrot), order);
216 for(i=0; i<nSH; i++){
217 for(j=0; j<nSH; j++) {
218 if(j==i)
219 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 1.0f, Mrot[i][j]);
220 else
221 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 0.0f, Mrot[i][j]);
222 }
223 }
224 free(Mrot);
225
226 /* Compare to the getSHrotMtx() Matlab function */
227 order = 4;
228 nSH = ORDER2NSH(order);
229 Mrot = (float**)malloc2d(nSH, nSH, sizeof(float));
230 yawPitchRoll2Rzyx(0.04f, 0.54f,-0.4f, 0, Rzyx);
231 getSHrotMtxReal(Rzyx, FLATTEN2D(Mrot), order);
232 double Mrot_ref[25][25] = {
233 {1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},{0,0.912317819470322,-0.334007492880439,-0.236886451652771,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
234 {0,0.408043822669133,0.790002010621868,0.457599237319041,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
235 {0,0.0342991990938353,-0.514135991653113,0.857022605902780,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
236 {0,0,0,0,0.773751979486127,-0.480511616313319,0.297436898769771,-0.164460121209763,-0.234308814625387,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
237 {0,0,0,0,0.320815885111266,0.584443217512645,-0.457030341925157,-0.339982347095703,-0.480664710153360,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
238 {0,0,0,0,0.323409465640717,0.558336000748573,0.436154765179890,0.626143845136656,0.0371501522262563,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
239 {0,0,0,0,0.365398067572425,-0.182693579159072,-0.703504421517165,0.441781344152855,0.378177314513551,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
240 {0,0,0,0,0.245510920021695,0.287086534852415,0.132306868781138,-0.519748017168846,0.754759962358177,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
241 {0,0,0,0,0,0,0,0,0,0.642754542747763,-0.587652464622319,0.146359326676735,-0.179940097166632,0.249957116297551,-0.161211805496773,-0.315061710316419,0,0,0,0,0,0,0,0,0},
242 {0,0,0,0,0,0,0,0,0,0.316547622267400,0.324276933833715,-0.489415761677808,0.525421745728824,-0.0811795764406443,-0.0642914639380568,-0.517998801533831,0,0,0,0,0,0,0,0,0},
243 {0,0,0,0,0,0,0,0,0,-0.0477608186606479,0.302122638638019,0.214473275742620,-0.433723919089070,-0.427443247772927,-0.611726955971008,-0.339717518973177,0,0,0,0,0,0,0,0,0},
244 {0,0,0,0,0,0,0,0,0,0.148935636035543,0.571302238306694,0.529863460253249,0.0476038953094580,0.594213419796629,0.0656256769672685,-0.104948528910382,0,0,0,0,0,0,0,0,0},
245 {0,0,0,0,0,0,0,0,0,0.311309233760352,0.304630835298635,-0.396153335826512,-0.667628966408715,-0.0103234397880398,0.454946318162605,0.0231945482299087,0,0,0,0,0,0,0,0,0},
246 {0,0,0,0,0,0,0,0,0,0.514785682894208,0.113244732089517,0.407883773582348,0.233719845299723,-0.593950310633879,0.241281704427283,0.300305444687571,0,0,0,0,0,0,0,0,0},
247 {0,0,0,0,0,0,0,0,0,0.316675769196523,0.161927142796105,-0.298312669792114,0.0285933354722383,0.205549150173188,-0.571110978701303,0.644414328446904,0,0,0,0,0,0,0,0,0},
248 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.526471642263643,-0.616929911516989,0.267922897453092,0.0235630456100945,0.0776050535864247,-0.190481327947399,0.295565129451190,-0.0753134473777231,-0.366811472459093},
249 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.234144273956354,0.0978665390875757,-0.545910447747527,0.175528558261790,-0.376101588123769,0.335795191612168,-0.141736252789070,-0.0455702308901721,-0.574798644029333},
250 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.0718436126062899,0.305262278899232,-0.0197737560173443,-0.298299395229287,0.646776790379034,0.111401675977437,0.0997398996043224,-0.463839920427382,-0.395542458465569},
251 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-0.155033529872408,-0.118946002867737,0.138228495430813,-0.0977208017941514,-0.285522105871139,-0.450196541284017,-0.600496309285322,-0.520682311298467,-0.131355606942160},
252 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0236933293789157,0.311297649179989,0.703254159219873,0.348811131545197,-0.261303521121084,0.391172954707122,0.0807830377413570,-0.219358047572331,-0.101769931423874},
253 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.146767948839247,0.439950893376704,0.0598087344890290,-0.520771343866458,-0.439502688322895,-0.362741803354952,0.407296904607327,0.0826968395396408,-0.112466610956744},
254 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.386795790652846,0.451176951621299,0.0223488932476933,0.463808781391941,0.287701399151563,-0.482347736946315,-0.226762742725175,0.241251512069808,-0.0784553883303562},
255 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.576800968786616,0.0555128465726625,0.144555412279657,-0.473213285269062,0.0597643274078365,0.343735767588532,-0.480720100388111,0.108090832343090,0.234286982126144},
256 {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.366598721881537,0.0733558553140817,-0.301930038675134,0.195400170636906,-0.0699710544219968,-0.0214401526687090,0.258994980191915,-0.617374325026823,0.526589247038282}};
257 for(i=0; i<nSH; i++)
258 for(j=0; j<nSH; j++)
259 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, (float)Mrot_ref[i][j], Mrot[i][j]);
260 free(Mrot);
261}
262
263void test__real2complexSHMtx(void){
264 int o, it, j, nSH, order;
265 float* Y_real_ref;
266 float_complex* Y_complex_ref, * tmp, *Y_complex_test;
267 float_complex** T_r2c;
268 float dir[2];
269 const float_complex calpha = cmplxf(1.0f, 0.0f), cbeta = cmplxf(0.0f, 0.0f);
270
271 /* Config */
272 const float acceptedTolerance = 0.0000001f;
273 int nTestOrders = 10;
274 int testOrders[10] = {1,2,3,4,5,6,7,8,9,10};
275 int nIter = 400;
276
277 /* Loop over orders */
278 for(o=0; o<nTestOrders; o++){
279 order = testOrders[o];
280 nSH = ORDER2NSH(order);
281 Y_real_ref = malloc1d(nSH * sizeof(float));
282 tmp = malloc1d(nSH * sizeof(float_complex));
283 Y_complex_ref = malloc1d(nSH * sizeof(float_complex));
284 Y_complex_test = malloc1d(nSH * sizeof(float_complex));
285 T_r2c = (float_complex**)malloc2d(nSH, nSH, sizeof(float_complex));
286
287 /* Loop over iterations */
288 for(it=0; it<nIter; it++) {
289 /* Random direction */
290 rand_m1_1(&dir[0] , 1);
291 rand_m1_1(&dir[1] , 1);
292 dir[0] *= SAF_PI;
293 dir[1] *= SAF_PI/2.0f;
294
295 /* Compute reference spherical harmonic weights */
296 getSHcomplex(order, (float*)dir, 1, Y_complex_ref);
297 getSHreal(order, (float*)dir, 1, Y_real_ref);
298
299 /* Convert to complex weights */
300 real2complexSHMtx(order, FLATTEN2D(T_r2c));
301 for(j=0; j<nSH; j++)
302 tmp[j] = cmplxf(Y_real_ref[j], 0.0f);
303 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasTrans, 1, nSH, nSH, &calpha,
304 tmp, nSH,
305 FLATTEN2D(T_r2c), nSH, &cbeta, /* Had to transpose it! */
306 Y_complex_test, nSH);
307
308 /* Should be equal to the reference */
309 for (j = 0; j < nSH; j++) {
310 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, crealf(Y_complex_ref[j]), crealf(Y_complex_test[j]));
311 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, cimagf(Y_complex_ref[j]), cimagf(Y_complex_test[j]));
312 }
313 }
314
315 /* Clean-up */
316 free(tmp);
317 free(Y_real_ref);
318 free(Y_complex_ref);
319 free(Y_complex_test);
320 free(T_r2c);
321 }
322}
323
324void test__complex2realSHMtx(void){
325 int o, it, j, nSH, order;
326 float* Y_real_ref;
327 float_complex* Y_complex_ref, *Y_real_test;
328 float_complex** T_c2r;
329 float dir[2];
330 const float_complex calpha = cmplxf(1.0f, 0.0f), cbeta = cmplxf(0.0f, 0.0f);
331
332 /* Config */
333 const float acceptedTolerance = 0.000001f;
334 int nTestOrders = 10;
335 int testOrders[10] = {1,2,3,4,5,6,7,8,9,10};
336 int nIter = 400;
337
338 /* Loop over orders */
339 for(o=0; o<nTestOrders; o++){
340 order = testOrders[o];
341 nSH = ORDER2NSH(order);
342 Y_real_ref = malloc1d(nSH * sizeof(float));
343 Y_complex_ref = malloc1d(nSH * sizeof(float_complex));
344 Y_real_test = malloc1d(nSH * sizeof(float_complex));
345 T_c2r = (float_complex**)malloc2d(nSH, nSH, sizeof(float_complex));
346
347 /* Loop over iterations */
348 for(it=0; it<nIter; it++) {
349 /* Random direction */
350 rand_m1_1(&dir[0] , 1);
351 rand_m1_1(&dir[1] , 1);
352 dir[0] *= SAF_PI;
353 dir[1] *= SAF_PI/2.0f;
354
355 /* Compute reference spherical harmonic weights */
356 getSHcomplex(order, (float*)dir, 1, Y_complex_ref);
357 getSHreal(order, (float*)dir, 1, Y_real_ref);
358
359 /* Convert to complex weights */
360 complex2realSHMtx(order, FLATTEN2D(T_c2r));
361 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasTrans, 1, nSH, nSH, &calpha,
362 Y_complex_ref, nSH,
363 FLATTEN2D(T_c2r), nSH, &cbeta, /* Had to transpose it! */
364 Y_real_test, nSH);
365
366 /* Should be equal to the reference */
367 for (j = 0; j < nSH; j++) {
368 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, Y_real_ref[j], crealf(Y_real_test[j]));
369 }
370 }
371
372 /* Clean-up */
373 free(Y_real_ref);
374 free(Y_complex_ref);
375 free(Y_real_test);
376 free(T_c2r);
377 }
378}
379
380void test__computeSectorCoeffsEP(void){
381 int i, j, numSec, order_sec, nSH_sec, nSH;
382 float* sec_dirs_deg;
383 float** sectorCoeffs;
384 float_complex*** A_xyz;
385
386 /* Config */
387 const float acceptedTolerance = 0.000001f;
388 const int order = 2;
389
390 /* Sector design and compute coefficients */
391 order_sec = order-1;
392 numSec = __Tdesign_nPoints_per_degree[2*order_sec-1];
393 sec_dirs_deg = (float*)__HANDLES_Tdesign_dirs_deg[2*order_sec-1];
394 nSH = ORDER2NSH(order);
395 nSH_sec = ORDER2NSH(order_sec);
396 A_xyz = (float_complex***)malloc3d(nSH, nSH_sec, 3, sizeof(float_complex));
397 computeVelCoeffsMtx(order_sec, FLATTEN3D(A_xyz));
398 sectorCoeffs = (float**)malloc2d((numSec*4),nSH,sizeof(float));
399 computeSectorCoeffsEP(order_sec, FLATTEN3D(A_xyz), SECTOR_PATTERN_PWD, sec_dirs_deg, numSec, FLATTEN2D(sectorCoeffs));
400
401 /* Check with Matlab reference */
402 double sectorCoeffs_ref[9][16]= {
403 {0.886226925452758,0.511663353973244,0.511663353973244,0.511663353973244,0.886226925452758,0.511663353973244,-0.511663353973244,-0.511663353973244,0.886226925452758,-0.511663353973244,0.511663353973244,-0.511663353973244,0.886226925452758,-0.511663353973244,-0.511663353973244,0.511663353973244},
404 {0.886226925452758,0,0.511663353973244,0,-0.886226925452758,0,0.511663353973244,0,0.886226925452758,0,0.511663353973244,0,-0.886226925452758,0,0.511663353973244,0},
405 {0.886226925452758,0,0,0.511663353973244,-0.886226925452758,0,0,0.511663353973244,-0.886226925452758,0,0,0.511663353973244,0.886226925452758,0,0,0.511663353973244},
406 {0.886226925452758,0.511663353973244,0,0,0.886226925452758,0.511663353973244,0,0,-0.886226925452758,0.511663353973244,0,0,-0.886226925452758,0.511663353973244,0,0},
407 {0,0.396332729760601,0.396332729760601,0,0,-0.396332729760601,0.396332729760601,0,0,0.396332729760601,-0.396332729760601,0,0,-0.396332729760601,-0.396332729760601,0},
408 {0,0,0.396332729760601,0.396332729760601,0,0,-0.396332729760601,-0.396332729760601,0,0,-0.396332729760601,0.396332729760601,0,0,0.396332729760601,-0.396332729760601},
409 {0,-0.228822808215942,-0.228822808215942,0.457645616431885,0,-0.228822808215942,0.228822808215942,-0.457645616431885,0,0.228822808215942,-0.228822808215942,-0.457645616431885,0,0.228822808215942,0.228822808215942,0.457645616431885},
410 {0,0.396332729760601,0,0.396332729760601,0,-0.396332729760601,0,0.396332729760601,0,-0.396332729760601,0,-0.396332729760601,0,0.396332729760601,0,-0.396332729760601},
411 {0,0.396332729760601,-0.396332729760601,0,0,0.396332729760601,0.396332729760601,0,0,-0.396332729760601,-0.396332729760601,0,0,-0.396332729760601,0.396332729760601,0}
412 };
413 for(i=0; i<9; i++)
414 for(j=0; j<16; j++)
415 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, (float)sectorCoeffs_ref[i][j], sectorCoeffs[j][i]);
416
417 free(sectorCoeffs);
418 free(A_xyz);
419}
420
421void test__checkCondNumberSHTReal(void){
422 int i, j, nDirs, order, nSH;
423 float* t_dirs_deg, *cond_N;
424 float** t_dirs_rad;
425
426 /* Config */
427 const float acceptedTolerance = 0.00001f;
428 int nTestOrders = 10;
429 int testOrders[10] = {1,2,3,4,5,6,7,8,9,10};
430
431 /* Loop over orders */
432 for(i=0; i<nTestOrders; i++) {
433 order = testOrders[i];
434 nSH = ORDER2NSH(order);
435
436 /* Pull an appropriate t-design */
437 t_dirs_deg = (float*)__HANDLES_Tdesign_dirs_deg[2 * order];
438 nDirs = __Tdesign_nPoints_per_degree[2 * order];
439 t_dirs_rad = (float **) malloc2d(nDirs, 2, sizeof(float));
440 for (j = 0; j < nDirs; j++) {
441 t_dirs_rad[j][0] = t_dirs_deg[j * 2] * SAF_PI / 180.0f;
442 t_dirs_rad[j][1] = SAF_PI / 2.0f - t_dirs_deg[j * 2 + 1] * SAF_PI /
443 180.0f; /* elevation->inclination */
444 }
445
446 /* Condition numbers for an appropriate t-design should be 1 */
447 cond_N = malloc1d((order+1)*sizeof(float));
448 checkCondNumberSHTReal(order, FLATTEN2D(t_dirs_rad), nDirs, NULL, cond_N);
449 for(j=0; j<order+1; j++)
450 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 1.0f, cond_N[j]);
451
452 /* Clean-up */
453 free(t_dirs_rad);
454 free(cond_N);
455 }
456}
457
458void test__calculateGridWeights(void){
459 int i, nDirs, order;
460 float* t_dirs_deg, *t_dirs_rad, *w;
461
462 /* Config */
463 const float acceptedTolerance = 0.00001f;
464 int testOrder = 3;
465
466 /* Pull an appropriate t-design */
467 t_dirs_deg = (float*)__HANDLES_Tdesign_dirs_deg[2 * testOrder -1];
468 nDirs = __Tdesign_nPoints_per_degree[2 * testOrder - 1];
469 t_dirs_rad = (float *) malloc1d(nDirs*2 *sizeof(float));
470 for(i=0; i<2*nDirs; i++)
471 t_dirs_rad[i] = DEG2RAD(t_dirs_deg[i]);
472 for(i=0; i<nDirs; i++)
473 t_dirs_rad[i*2+1] = SAF_PI/2.0f - t_dirs_rad[i*2+1];
474 w = malloc1d(nDirs * sizeof(float));
475 order = calculateGridWeights(t_dirs_rad,nDirs,-1,w);
476
477 TEST_ASSERT_EQUAL(testOrder, order);
478 for (i=0; i<nDirs; i++)
479 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance,FOURPI/nDirs, w[i]);
480
481 /* clean-up */
482 free(t_dirs_rad);
483 free(w);
484}
485
486void test__sphMUSIC(void){
487 int i, j, k, nGrid, nSH, nSrcs, srcInd_1, srcInd_2;
488 float test_dirs_deg[2][2];
489 float* grid_dirs_deg;
490 float** Y_src, **src_sigs, **src_sigs_sh, **Cx, **V, **Vn;
491 float_complex** Vn_cmplx;
492 void* hMUSIC;
493
494 /* config */
495 const int order = 3;
496 const int lsig = 48000;
497
498 /* define scanning grid directions */
499 nGrid = 240;
500 grid_dirs_deg = (float*)__Tdesign_degree_21_dirs_deg;
501
502 /* test scenario and signals */
503 nSrcs = 2;
504 srcInd_1 = 139;
505 srcInd_2 = 204;
506 test_dirs_deg[0][0] = grid_dirs_deg[srcInd_1*2];
507 test_dirs_deg[0][1] = grid_dirs_deg[srcInd_1*2+1];
508 test_dirs_deg[1][0] = grid_dirs_deg[srcInd_2*2];
509 test_dirs_deg[1][1] = grid_dirs_deg[srcInd_2*2+1];
510 nSH = ORDER2NSH(order);
511 Y_src = (float**)malloc2d(nSH, nSrcs, sizeof(float));
512 getRSH(order, (float*)test_dirs_deg, nSrcs, FLATTEN2D(Y_src));
513 src_sigs = (float**)malloc2d(nSrcs, lsig, sizeof(float));
514 rand_m1_1(FLATTEN2D(src_sigs), nSrcs*lsig); /* uncorrelated noise sources */
515
516 /* encode to SH and compute spatial covariance matrix */
517 src_sigs_sh = (float**)malloc2d(nSH, lsig, sizeof(float));
518 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nSH, lsig, nSrcs, 1.0f,
519 FLATTEN2D(Y_src), nSrcs,
520 FLATTEN2D(src_sigs), lsig, 0.0f,
521 FLATTEN2D(src_sigs_sh), lsig);
522 Cx = (float**)malloc2d(nSH, nSH, sizeof(float));
523 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, nSH, nSH, lsig, 1.0f,
524 FLATTEN2D(src_sigs_sh), lsig,
525 FLATTEN2D(src_sigs_sh), lsig, 0.0f,
526 FLATTEN2D(Cx), nSH);
527
528 /* Eigenvalue decomposition and truncation of eigen vectors to obtain
529 * noise subspace (based on source number) */
530 V = (float**)malloc2d(nSH, nSH, sizeof(float));
531 utility_sseig(NULL, FLATTEN2D(Cx), nSH, 1, FLATTEN2D(V), NULL, NULL);
532 Vn = (float**)malloc2d(nSH, (nSH-nSrcs), sizeof(float)); /* noise subspace */
533 for(i=0; i<nSH; i++)
534 for(j=0, k=nSrcs; j<nSH-nSrcs; j++, k++)
535 Vn[i][j] = V[i][k];
536 Vn_cmplx = (float_complex**)malloc2d(nSH, (nSH-nSrcs), sizeof(float_complex)); /* noise subspace (complex) */
537 for(i=0; i<nSH; i++)
538 for(j=0; j<nSH-nSrcs; j++)
539 Vn_cmplx[i][j] = cmplxf(Vn[i][j], 0.0f);
540
541 /* compute sphMUSIC, returning "peak-find" indices */
542 int inds[2];
543 sphMUSIC_create(&hMUSIC, order, grid_dirs_deg, nGrid);
544 sphMUSIC_compute(hMUSIC, FLATTEN2D(Vn_cmplx), nSrcs, NULL, (int*)inds);
545
546 /* Assert that the true source indices were found (note that the order can flip) */
547 TEST_ASSERT_TRUE(inds[0] == srcInd_1 || inds[0] == srcInd_2);
548 TEST_ASSERT_TRUE(inds[1] == srcInd_1 || inds[1] == srcInd_2);
549
550 /* clean-up */
551 sphMUSIC_destroy(&hMUSIC);
552 free(Y_src);
553 free(src_sigs);
554 free(src_sigs_sh);
555 free(Cx);
556 free(V);
557 free(Vn);
558 free(Vn_cmplx);
559}
560
561void test__sphPWD(void){
562 int nGrid, nSH, nSrcs, srcInd_1, srcInd_2;
563 float test_dirs_deg[2][2];
564 float* grid_dirs_deg;
565 float** Y_src, **src_sigs, **src_sigs_sh, **Cx;
566 float_complex** Cx_cmplx;
567 void* hPWD;
568
569 /* config */
570 const int order = 3;
571 const int lsig = 48000;
572
573 /* define scanning grid directions */
574 nGrid = 240;
575 grid_dirs_deg = (float*)__Tdesign_degree_21_dirs_deg;
576
577 /* test scenario and signals */
578 nSrcs = 2;
579 srcInd_1 = 139;
580 srcInd_2 = 204;
581 test_dirs_deg[0][0] = grid_dirs_deg[srcInd_1*2];
582 test_dirs_deg[0][1] = grid_dirs_deg[srcInd_1*2+1];
583 test_dirs_deg[1][0] = grid_dirs_deg[srcInd_2*2];
584 test_dirs_deg[1][1] = grid_dirs_deg[srcInd_2*2+1];
585 nSH = ORDER2NSH(order);
586 Y_src = (float**)malloc2d(nSH, nSrcs, sizeof(float));
587 getRSH(order, (float*)test_dirs_deg, nSrcs, FLATTEN2D(Y_src));
588 src_sigs = (float**)malloc2d(nSrcs, lsig, sizeof(float));
589 rand_m1_1(FLATTEN2D(src_sigs), nSrcs*lsig); /* uncorrelated noise sources */
590
591 /* encode to SH and compute spatial covariance matrix */
592 src_sigs_sh = (float**)malloc2d(nSH, lsig, sizeof(float));
593 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nSH, lsig, nSrcs, 1.0f,
594 FLATTEN2D(Y_src), nSrcs,
595 FLATTEN2D(src_sigs), lsig, 0.0f,
596 FLATTEN2D(src_sigs_sh), lsig);
597 Cx = (float**)malloc2d(nSH, nSH, sizeof(float));
598 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, nSH, nSH, lsig, 1.0f,
599 FLATTEN2D(src_sigs_sh), lsig,
600 FLATTEN2D(src_sigs_sh), lsig, 0.0f,
601 FLATTEN2D(Cx), nSH);
602 Cx_cmplx = (float_complex**)calloc2d(nSH, nSH, sizeof(float_complex));
603 cblas_scopy(nSH*nSH, FLATTEN2D(Cx), 1, (float*)FLATTEN2D(Cx_cmplx), 2);
604
605 /* compute sphMUSIC, returning "peak-find" indices */
606 int inds[2];
607 sphPWD_create(&hPWD, order, grid_dirs_deg, nGrid);
608 sphPWD_compute(hPWD, FLATTEN2D(Cx_cmplx), nSrcs, NULL, (int*)inds);
609
610 /* Assert that the true source indices were found (note that the order can flip) */
611 TEST_ASSERT_TRUE(inds[0] == srcInd_1 || inds[0] == srcInd_2);
612 TEST_ASSERT_TRUE(inds[1] == srcInd_1 || inds[1] == srcInd_2);
613
614 /* clean-up */
615 sphPWD_destroy(&hPWD);
616 free(Y_src);
617 free(src_sigs);
618 free(src_sigs_sh);
619 free(Cx);
620 free(Cx_cmplx);
621}
622
623void test__sphESPRIT(void){
624 int i,j,nSH, nSrcs;
625 void* hESPRIT;
626 float test_dirs_deg[2][2], estdirs_deg[2][2];
627 float** Y_src, **src_sigs, **src_sigs_sh, **tmpCx;
628 float_complex** Cx, **C_Cx, **T_r2c, **Cx_R, **U, **Us;
629 const float_complex calpha = cmplxf(1.0f, 0.0f), cbeta = cmplxf(0.0f, 0.0f);
630
631 /* config */
632 const float acceptedTolerance = 0.01f; /* degrees */
633 const int order = 3;
634 const int lsig = 48000;
635
636 /* test scenario and signals */
637 nSrcs = 2;
638 test_dirs_deg[0][0] = -90.0f;
639 test_dirs_deg[0][1] = 10.0f;
640 test_dirs_deg[1][0] = 20.0f;
641 test_dirs_deg[1][1] = -40.0f;
642 nSH = ORDER2NSH(order);
643 Y_src = (float**)malloc2d(nSH, nSrcs, sizeof(float));
644 getRSH(order, (float*)test_dirs_deg, nSrcs, FLATTEN2D(Y_src));
645 src_sigs = (float**)malloc2d(nSrcs, lsig, sizeof(float));
646 rand_m1_1(FLATTEN2D(src_sigs), nSrcs*lsig); /* uncorrelated noise sources */
647
648 /* encode to SH and compute spatial covariance matrix */
649 src_sigs_sh = (float**)malloc2d(nSH, lsig, sizeof(float));
650 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nSH, lsig, nSrcs, 1.0f,
651 FLATTEN2D(Y_src), nSrcs,
652 FLATTEN2D(src_sigs), lsig, 0.0f,
653 FLATTEN2D(src_sigs_sh), lsig);
654 tmpCx = (float**)malloc2d(nSH, nSH, sizeof(float));
655 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, nSH, nSH, lsig, 1.0f,
656 FLATTEN2D(src_sigs_sh), lsig,
657 FLATTEN2D(src_sigs_sh), lsig, 0.0f,
658 FLATTEN2D(tmpCx), nSH);
659 Cx = (float_complex**)malloc2d(nSH, nSH, sizeof(float_complex));
660 for(i=0; i<nSH; i++)
661 for(j=0; j<nSH; j++)
662 Cx[i][j] = cmplxf(tmpCx[i][j], 0.0f); /* real->complex data-type */
663
664 /* Convert to complex basis */
665 T_r2c = (float_complex**)malloc2d(nSH, nSH, sizeof(float_complex));
666 real2complexSHMtx(order, FLATTEN2D(T_r2c));
667 for(i=0; i<nSH; i++)
668 for(j=0; j<nSH; j++)
669 T_r2c[i][j] = conjf(T_r2c[i][j]);
670 Cx_R = (float_complex**)malloc2d(nSH, nSH, sizeof(float_complex));
671 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, nSH, nSH, nSH, &calpha,
672 FLATTEN2D(Cx), nSH,
673 FLATTEN2D(T_r2c), nSH, &cbeta,
674 FLATTEN2D(Cx_R), nSH);
675 C_Cx = (float_complex**)malloc2d(nSH, nSH, sizeof(float_complex));
676 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nSH, nSH, nSH, &calpha,
677 FLATTEN2D(T_r2c), nSH,
678 FLATTEN2D(Cx_R), nSH, &cbeta,
679 FLATTEN2D(C_Cx), nSH);
680
681 /* Eigenvalue decomposition and truncation of eigen vectors to obtain
682 * signal subspace (based on source number) */
683 U = (float_complex**)malloc2d(nSH, nSH, sizeof(float_complex));
684 utility_cseig(NULL, FLATTEN2D(C_Cx), nSH, 1, FLATTEN2D(U), NULL, NULL);
685 Us = (float_complex**)malloc2d(nSH, nSrcs, sizeof(float_complex)); /* signal subspace */
686 for(i=0; i<nSH; i++)
687 for(j=0; j<nSrcs; j++)
688 Us[i][j] = U[i][j];
689
690 /* use sphESPRIT to estimate source directions... */
691 sphESPRIT_create(&hESPRIT, order);
692 sphESPRIT_estimateDirs(hESPRIT, FLATTEN2D(Us), nSrcs, (float*)estdirs_deg);
693 for(i=0; i<nSrcs; i++){
694 estdirs_deg[i][0]*=180.0f/SAF_PI; /* rad->deg */
695 estdirs_deg[i][1]*=180.0f/SAF_PI;
696 }
697
698 /* Assert that the true source directions were found (note that the order can flip) */
699 TEST_ASSERT_TRUE(fabsf(estdirs_deg[0][0]-test_dirs_deg[0][0])<acceptedTolerance ||
700 fabsf(estdirs_deg[1][0]-test_dirs_deg[0][0])<acceptedTolerance);
701 TEST_ASSERT_TRUE(fabsf(estdirs_deg[0][1]-test_dirs_deg[0][1])<acceptedTolerance ||
702 fabsf(estdirs_deg[1][1]-test_dirs_deg[0][1])<acceptedTolerance);
703 TEST_ASSERT_TRUE(fabsf(estdirs_deg[0][0]-test_dirs_deg[1][0])<acceptedTolerance ||
704 fabsf(estdirs_deg[1][0]-test_dirs_deg[1][0])<acceptedTolerance);
705 TEST_ASSERT_TRUE(fabsf(estdirs_deg[0][1]-test_dirs_deg[1][1])<acceptedTolerance ||
706 fabsf(estdirs_deg[1][1]-test_dirs_deg[1][1])<acceptedTolerance);
707
708 /* clean-up */
709 sphESPRIT_destroy(&hESPRIT);
710 free(Y_src);
711 free(src_sigs);
712 free(src_sigs_sh);
713 free(tmpCx);
714 free(Cx);
715 free(C_Cx);
716 free(T_r2c);
717 free(Cx_R);
718 free(U);
719 free(Us);
720}
721
722void test__sphModalCoeffs(void){
723 int i, j;
724 float* freqVector;
725 double* kr;
726 double_complex** b_N_dipole, **b_N_card, **b_N_omni, **b_N_omni_test;
727
728 /* Config */
729 const double acceptedTolerance = 0.000001f;
730 const int order = 4;
731 const int N = 16;
732 const float fs = 48000;
733 const double radius = 0.04;
734 const double c = 343.0;
735
736 /* prep */
737 freqVector = malloc1d((N/2+1)*sizeof(float));
738 getUniformFreqVector(N, fs, freqVector);
739 kr = malloc1d((N/2+1)*sizeof(double));
740 for(i=0; i<N/2+1; i++)
741 kr[i] = 2.0*SAF_PId* (double)freqVector[i] * radius/c;
742 b_N_dipole = (double_complex**)malloc2d((N/2+1), (order+1), sizeof(double_complex));
743 b_N_card = (double_complex**)malloc2d((N/2+1), (order+1), sizeof(double_complex));
744 b_N_omni = (double_complex**)malloc2d((N/2+1), (order+1), sizeof(double_complex));
745 b_N_omni_test = (double_complex**)malloc2d((N/2+1), (order+1), sizeof(double_complex));
746
747 /* Compute modal coefficients */
748 sphModalCoeffs(order, kr, (N/2+1), ARRAY_CONSTRUCTION_OPEN_DIRECTIONAL, 0.0, FLATTEN2D(b_N_dipole));
749 sphModalCoeffs(order, kr, (N/2+1), ARRAY_CONSTRUCTION_OPEN_DIRECTIONAL, 0.5, FLATTEN2D(b_N_card));
750 sphModalCoeffs(order, kr, (N/2+1), ARRAY_CONSTRUCTION_OPEN_DIRECTIONAL, 1.0, FLATTEN2D(b_N_omni));
751 sphModalCoeffs(order, kr, (N/2+1), ARRAY_CONSTRUCTION_OPEN, 666.0 /* not used */, FLATTEN2D(b_N_omni_test));
752
753 /* Check that "open directional", with "dirCoeff=1" is identical to just "open" */
754 for(i=0; i<N/2+1; i++){
755 for(j=0; j< order+1; j++){
756 TEST_ASSERT_TRUE( fabs(creal(b_N_omni[i][j]) - creal(b_N_omni_test[i][j])) <= acceptedTolerance );
757 TEST_ASSERT_TRUE( fabs(cimag(b_N_omni[i][j]) - cimag(b_N_omni_test[i][j])) <= acceptedTolerance );
758 }
759 }
760
761 /* clean-up */
762 free(b_N_dipole);
763 free(b_N_card);
764 free(b_N_omni);
765 free(b_N_omni_test);
766}
767
768void test__arraySHTmatrices(void){
769 int nGrid;
770 float* grid_dirs_deg;
771 float_complex* H_array, *H_sht;
772
773 /* Configuration */
774 const int order = 2;
775 const int nMics = 13;
776 const int nBins = 129;
777 const float amp_thresh_dB = 15.0f;
778
779 /* Prep */
780 nGrid = __Tdesign_degree_30_nPoints;
781 grid_dirs_deg = (float*)__Tdesign_degree_30_dirs_deg;
782 H_array = malloc1d(nBins*nMics*nGrid*sizeof(float_complex));
783 rand_m1_1((float*)H_array, /*re+im*/2*nBins*nMics*nGrid);
784 H_sht = malloc1d(nBins*ORDER2NSH(order)*nMics*sizeof(float_complex));
785
786 /* (Tikhonov) regularised LS solution */
787 arraySHTmatrices(ARRAY_SHT_REG_LS, order, amp_thresh_dB, H_array, grid_dirs_deg, nBins, nMics, nGrid, NULL, H_sht);
788
789 /* (Tikhonov) regularised LS solution in the spherical harmonic domain */
790 arraySHTmatrices(ARRAY_SHT_REG_LSHD, order, amp_thresh_dB, H_array, grid_dirs_deg, nBins, nMics, nGrid, NULL, H_sht);
791
792 /* clean-up */
793 free(H_array);
794 free(H_sht);
795}
getRSH
void getRSH(int N, float *dirs_deg, int nDirs, float *Y)
Computes real-valued spherical harmonics [1] for each given direction on the unit sphere.
Definition saf_hoa.c:119
sphMUSIC_destroy
void sphMUSIC_destroy(void **const phMUSIC)
Destroys an instance of the spherical harmonic domain MUSIC implementation.
Definition saf_sh.c:1333
ORDER2NSH
#define ORDER2NSH(order)
Converts spherical harmonic order to number of spherical harmonic components i.e: (order+1)^2.
Definition saf_sh.h:51
sphPWD_create
void sphPWD_create(void **const phPWD, int order, float *grid_dirs_deg, int nDirs)
Creates an instance of a spherical harmonic domain implementation of the steer-response power (SRP) a...
Definition saf_sh.c:1155
sphMUSIC_compute
void sphMUSIC_compute(void *const hMUSIC, float_complex *Vn, int nSrcs, float *P_music, int *peak_inds)
Computes a pseudo-spectrum based on the MUSIC algorithm in the spherical harmonic domain; optionally ...
Definition saf_sh.c:1356
sphESPRIT_create
void sphESPRIT_create(void **const phESPRIT, int order)
Creates an instance of the spherical harmonic domain ESPRIT-based direction of arrival estimator.
Definition saf_sh.c:1421
computeSectorCoeffsEP
float computeSectorCoeffsEP(int orderSec, float_complex *A_xyz, SECTOR_PATTERNS pattern, float *sec_dirs_deg, int nSecDirs, float *sectorCoeffs)
Computes beamforming matrices (sector coefficients) which, when applied to input SH signals,...
Definition saf_sh.c:704
calculateGridWeights
int calculateGridWeights(float *dirs_rad, int nDirs, int order, float *w)
Computes approximation of quadrature/integration weights for a given grid.
Definition saf_sh.c:1066
getSHcomplex
void getSHcomplex(int order, float *dirs_rad, int nDirs, float_complex *Y)
Computes complex-valued spherical harmonics [1] for each given direction on the unit sphere.
Definition saf_sh.c:440
sphMUSIC_create
void sphMUSIC_create(void **const phMUSIC, int order, float *grid_dirs_deg, int nDirs)
Creates an instance of the spherical harmonic domain MUSIC implementation, which may be used for comp...
Definition saf_sh.c:1285
sphESPRIT_estimateDirs
void sphESPRIT_estimateDirs(void *const hESPRIT, float_complex *Us, int K, float *src_dirs_rad)
Estimates the direction-of-arrival (DoA) based on the ESPRIT-based estimator, in the spherical harmon...
Definition saf_sh.c:1546
sphModalCoeffs
void sphModalCoeffs(int order, double *kr, int nBands, ARRAY_CONSTRUCTION_TYPES arrayType, double dirCoeff, double_complex *b_N)
Calculates the modal coefficients for open/rigid spherical arrays.
Definition saf_sh.c:2370
sphPWD_destroy
void sphPWD_destroy(void **const phPWD)
Destroys an instance of the spherical harmonic domain PWD implementation.
Definition saf_sh.c:1201
complex2realSHMtx
void complex2realSHMtx(int order, float_complex *T_c2r)
Computes a complex to real spherical harmonic transform matrix.
Definition saf_sh.c:491
sphPWD_compute
void sphPWD_compute(void *const hPWD, float_complex *Cx, int nSrcs, float *P_map, int *peak_inds)
Computes a power-map based on determining the energy of hyper-cardioid beamformers; optionally,...
Definition saf_sh.c:1222
arraySHTmatrices
void arraySHTmatrices(ARRAY_SHT_OPTIONS method, int order, float amp_thresh_dB, float_complex *H_array, float *grid_dirs_deg, int nBins, int nMics, int nGrid, float *w_grid, float_complex *H_sht)
Computes matrices required to transform array signals into spherical harmonic signals (frequency-doma...
Definition saf_sh.c:1978
sphESPRIT_destroy
void sphESPRIT_destroy(void **const phESPRIT)
Destroys an instance of the spherical harmonic domain ESPRIT-based direction of arrival estimator.
Definition saf_sh.c:1496
getSHrotMtxReal
void getSHrotMtxReal(float Rxyz[3][3], float *RotMtx, int L)
Generates a real-valued spherical harmonic rotation matrix [1] based on a 3x3 rotation matrix (see qu...
Definition saf_sh.c:586
getSHreal_recur
void getSHreal_recur(int N, float *dirs_rad, int nDirs, float *Y)
Computes real-valued spherical harmonics [1] for each given direction on the unit sphere.
Definition saf_sh.c:254
computeVelCoeffsMtx
void computeVelCoeffsMtx(int sectorOrder, float_complex *A_xyz)
Computes the matrices which generate the coefficients of a beampattern of order (sectorOrder+1) that ...
Definition saf_sh.c:672
getSHreal
void getSHreal(int order, float *dirs_rad, int nDirs, float *Y)
Computes real-valued spherical harmonics [1] for each given direction on the unit sphere.
Definition saf_sh.c:191
checkCondNumberSHTReal
void checkCondNumberSHTReal(int order, float *dirs_rad, int nDirs, float *w, float *cond_N)
Computes the condition numbers for a least-squares SHT.
Definition saf_sh.c:991
real2complexSHMtx
void real2complexSHMtx(int order, float_complex *T_r2c)
Computes a real to complex spherical harmonic transform matrix.
Definition saf_sh.c:523
ARRAY_CONSTRUCTION_OPEN_DIRECTIONAL
@ ARRAY_CONSTRUCTION_OPEN_DIRECTIONAL
Open array, directional sensors.
Definition saf_sh.h:67
ARRAY_CONSTRUCTION_OPEN
@ ARRAY_CONSTRUCTION_OPEN
Open array, omni-directional sensors.
Definition saf_sh.h:65
SECTOR_PATTERN_PWD
@ SECTOR_PATTERN_PWD
Plane-wave decomposition/Hyper-cardioid.
Definition saf_sh.h:82
ARRAY_SHT_REG_LS
@ ARRAY_SHT_REG_LS
Regularised least-squares (LS)
Definition saf_sh.h:100
ARRAY_SHT_REG_LSHD
@ ARRAY_SHT_REG_LSHD
Regularised least-squares (LS) in the SH domain (similar as in [1])
Definition saf_sh.h:101
SAF_PI
#define SAF_PI
pi constant (single precision)
Definition saf_utilities.h:70
__HANDLES_Tdesign_dirs_deg
const float * __HANDLES_Tdesign_dirs_deg[21]
minimum T-design HANDLES (up to degree 21 only).
Definition saf_utility_loudspeaker_presets.c:29910
__Tdesign_degree_30_nPoints
const int __Tdesign_degree_30_nPoints
Number of directions in a minimum Tdesign of degree: 30.
Definition saf_utility_loudspeaker_presets.c:2542
DEG2RAD
#define DEG2RAD(x)
Converts degrees to radians.
Definition saf_utilities.h:99
__Tdesign_degree_21_dirs_deg
const float __Tdesign_degree_21_dirs_deg[240][2]
Directions [azimuth, Elevation] in degrees, for minimum Tdesign degree: 21.
Definition saf_utility_loudspeaker_presets.c:2299
SAF_PId
#define SAF_PId
pi constant (double precision)
Definition saf_utilities.h:73
FOURPI
#define FOURPI
4pi (single precision)
Definition saf_utilities.h:92
__Tdesign_degree_30_dirs_deg
const float __Tdesign_degree_30_dirs_deg[480][2]
Directions [azimuth, Elevation] in degrees, for minimum Tdesign degree: 30.
Definition saf_utility_loudspeaker_presets.c:2543
utility_cseig
void utility_cseig(void *const hWork, const float_complex *A, const int dim, int sortDecFLAG, float_complex *V, float_complex *D, float *eig)
Eigenvalue decomposition of a SYMMETRIC/HERMITION matrix: single precision complex,...
Definition saf_utility_veclib.c:2096
getUniformFreqVector
void getUniformFreqVector(int fftSize, float fs, float *freqVector)
Calculates the frequencies (in Hz) of uniformly spaced bins, for a given FFT size and sampling rate.
Definition saf_utility_fft.c:146
utility_sseig
void utility_sseig(void *const hWork, const float *A, const int dim, int sortDecFLAG, float *V, float *D, float *eig)
Eigenvalue decomposition of a SYMMETRIC matrix: single precision, i.e.
Definition saf_utility_veclib.c:1955
__Tdesign_nPoints_per_degree
const int __Tdesign_nPoints_per_degree[21]
Number of points in each t-design (up to degree 21 only).
Definition saf_utility_loudspeaker_presets.c:29906
yawPitchRoll2Rzyx
void yawPitchRoll2Rzyx(float yaw, float pitch, float roll, int rollPitchYawFLAG, float R[3][3])
Constructs a 3x3 rotation matrix from the Euler angles, using the yaw-pitch-roll (zyx) convention.
Definition saf_utility_geometry.c:258
utility_cvsmul
void utility_cvsmul(float_complex *a, const float_complex *s, const int len, float_complex *c)
Single-precision, complex, multiplies each element in vector 'a' with a scalar 's',...
Definition saf_utility_veclib.c:1272
SQRT4PI
#define SQRT4PI
sqrt(4pi) (single precision)
Definition saf_utilities.h:89
utility_svsmul
void utility_svsmul(float *a, const float *s, const int len, float *c)
Single-precision, multiplies each element in vector 'a' with a scalar 's', i.e.
Definition saf_utility_veclib.c:1244
rand_m1_1
void rand_m1_1(float *vector, int length)
Generates random numbers between -1 and 1 and stores them in the input vector.
Definition saf_utility_misc.c:228
malloc2d
void ** malloc2d(size_t dim1, size_t dim2, size_t data_size)
2-D malloc (contiguously allocated, so use free() as usual to deallocate)
Definition md_malloc.c:89
malloc1d
void * malloc1d(size_t dim1_data_size)
1-D malloc (same as malloc, but with error checking)
Definition md_malloc.c:59
malloc3d
void *** malloc3d(size_t dim1, size_t dim2, size_t dim3, size_t data_size)
3-D malloc (contiguously allocated, so use free() as usual to deallocate)
Definition md_malloc.c:151
calloc2d
void ** calloc2d(size_t dim1, size_t dim2, size_t data_size)
2-D calloc (contiguously allocated, so use free() as usual to deallocate)
Definition md_malloc.c:102
FLATTEN2D
#define FLATTEN2D(A)
Use this macro when passing a 2-D dynamic multi-dimensional array to memset, memcpy or any other func...
Definition md_malloc.h:65
FLATTEN3D
#define FLATTEN3D(A)
Use this macro when passing a 3-D dynamic multi-dimensional array to memset, memcpy or any other func...
Definition md_malloc.h:72
saf_test.h
Unit test program for the Spatial_Audio_Framework.
test__real2complexSHMtx
void test__real2complexSHMtx(void)
Testing the real to complex spherical harmonic conversion, using getSHcomplex() as the reference.
Definition test__sh_module.c:263
test__calculateGridWeights
void test__calculateGridWeights(void)
Test grid weight approximation.
Definition test__sh_module.c:458
test__sphPWD
void test__sphPWD(void)
Testing the DoA estimation performance of sphPWD()
Definition test__sh_module.c:561
test__getSHreal_recur
void test__getSHreal_recur(void)
Testing that the getSHreal_recur() function is somewhat numerically identical to the full-fat getSHre...
Definition test__sh_module.c:84
test__getSHreal
void test__getSHreal(void)
Testing the orthogonality of the getSHreal() function.
Definition test__sh_module.c:27
test__arraySHTmatrices
void test__arraySHTmatrices(void)
Testing the arraySHTmatrices() function.
Definition test__sh_module.c:768
test__getSHreal_part
void test__getSHreal_part(void)
Testing that the getSHreal_part() function is somewhat numerically identical to the full-fat getSHrea...
Definition test__sh_module.c:111
test__sphESPRIT
void test__sphESPRIT(void)
Testing the DoA estimation performance of sphESPRIT()
Definition test__sh_module.c:623
test__getSHrotMtxReal
void test__getSHrotMtxReal(void)
Testing the spherical harmonic rotation matrix function getSHrotMtxReal()
Definition test__sh_module.c:202
test__sphModalCoeffs
void test__sphModalCoeffs(void)
Testing the sphModalCoeffs() function.
Definition test__sh_module.c:722
test__complex2realSHMtx
void test__complex2realSHMtx(void)
Testing the complex to real spherical harmonic conversion, using getSHreal() as the reference.
Definition test__sh_module.c:324
test__computeSectorCoeffsEP
void test__computeSectorCoeffsEP(void)
Testing the computeSectorCoeffsEP() and computeVelCoeffsMtx() functions and comparing their output to...
Definition test__sh_module.c:380
test__checkCondNumberSHTReal
void test__checkCondNumberSHTReal(void)
Testing that for T-designs, the condition numbers are all equal to 1.
Definition test__sh_module.c:421
test__getSHcomplex
void test__getSHcomplex(void)
Testing the orthogonality of the getSHcomplex() function.
Definition test__sh_module.c:143
test__sphMUSIC
void test__sphMUSIC(void)
Testing the DoA estimation performance of sphMUSIC()
Definition test__sh_module.c:486