SAF
Loading...
Searching...
No Matches
test__utilities_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"
26
27void test__cylindricalBesselFunctions(void){ // TODO: may as well check the derivatives too...
28 int i;
29 double J_n[10], Y_n[10];
30
31 /* Config */
32 const float acceptedTolerance = 0.00001f;
33 int testOrder = 7; /* note, REF values hardcoded for order 7 */
34 double z[10] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}; /* note, REF values hardcoded for these values */
35
36 /* Reference values computed in MATLAB with:
37 * J_n = besselj(N, z);
38 * y_n = bessely(N, z);
39 */
40 double J_nREF[10] = {0.0, 1.50232581743681e-06, 0.000174944074868274, 0.00254729445180469, 0.0151760694220584, 0.0533764101558907, 0.129586651841481, 0.233583569505696, 0.320589077979826, 0.327460879242453};
41 double Y_nREF[10] = {0.0, -30588.9570521240, -271.548025367994, -19.8399354089864, -3.70622393164077, -1.26289883576932, -0.656590825719075, -0.405371018606768, -0.200063904600409, 0.0172445799076681};
42
43 /* test bessel_Jn */
44 bessel_Jn(testOrder, z, 10, J_n, NULL);
45 for(i=0; i<10; i++)
46 TEST_ASSERT_TRUE(fabs(J_n[i]-J_nREF[i])<acceptedTolerance);
47
48 /* test bessel_Yn */
49 bessel_Yn(testOrder, z, 10, Y_n, NULL);
50 for(i=0; i<10; i++)
51 TEST_ASSERT_TRUE(fabs(Y_n[i]-Y_nREF[i])<acceptedTolerance);
52}
53
54void test__sphericalBesselFunctions(void){ // TODO: may as well check the derivatives too...
55 int i, successFlag;
56 double j_n[10], i_n[10], y_n[10], k_n[10];
57
58 /* Config */
59 const float acceptedTolerance = 0.00001f;
60 int testOrder = 7; /* note, REF values hardcoded for order 7 */
61 double z[10] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0}; /* note, REF values hardcoded for these values */
62
63 /* Reference values computed in MATLAB with:
64 * j_n = sqrt(pi./(2*z)).*besselj(N+0.5, z);
65 * i_n = sqrt(pi./(2*z)).*besseli(N+0.5, z);
66 * y_n = sqrt(pi./(2*z)).*bessely(N+0.5, z);
67 * k_n = sqrt(pi./(2*z)).*besselk(N+0.5, z);
68 */
69 double j_nREF[10] = {0.0, 4.79013419873948e-07, 5.60965570334894e-05, 0.000824843253217635, 0.00498650846172602, 0.0179027781779895, 0.0447223808293482, 0.0839226228445072, 0.122272711565833, 0.137946585027486};
70 double i_nREF[10] = {0.0, 5.08036087257580e-07, 7.09794452304064e-05, 0.00140087680258227, 0.0127983365433790, 0.0783315436379810, 0.377879458299915, 1.56419501808402, 5.83626393050750, 20.2384754394417};
71 double y_nREF[10] = {0.0, -140452.852366906, -617.054329642527, -29.4761692244538, -3.98778927238432, -1.02739463881260, -0.425887203702750, -0.237025274765842, -0.132622247946352, -0.0402143438632017};
72 double k_nREF[10] = {0.0, 204287.522076393, 712.406907885478, 23.1153112578315, 1.80293583642309, 0.222213613092395, 0.0360276414091966, 0.00698538879470478, 0.00153285534574965, 0.000367847412220325};
73
74 /* test bessel_jn */
75 successFlag = bessel_jn(testOrder, z, 10, j_n, NULL);
76 TEST_ASSERT_TRUE(successFlag);
77 for(i=0; i<10; i++)
78 TEST_ASSERT_TRUE(fabs(j_n[i]-j_nREF[i])<acceptedTolerance);
79
80 /* test bessel_in */
81 successFlag = bessel_in(testOrder, z, 10, i_n, NULL);
82 TEST_ASSERT_TRUE(successFlag);
83 for(i=0; i<10; i++)
84 TEST_ASSERT_TRUE(fabs(i_n[i]-i_nREF[i])<acceptedTolerance);
85
86 /* test bessel_yn */
87 successFlag = bessel_yn(testOrder, z, 10, y_n, NULL);
88 TEST_ASSERT_TRUE(successFlag);
89 for(i=0; i<10; i++)
90 TEST_ASSERT_TRUE(fabs(y_n[i]-y_nREF[i])<acceptedTolerance);
91
92 /* test bessel_kn */
93 successFlag = bessel_kn(testOrder, z, 10, k_n, NULL);
94 TEST_ASSERT_TRUE(successFlag);
95 for(i=0; i<10; i++)
96 TEST_ASSERT_TRUE(fabs(k_n[i]-k_nREF[i])<acceptedTolerance);
97}
98
99void test__cart2sph(void){
100 const float acceptedTolerance = 0.00001f;
101 float cordCar [100][3];
102 float cordSph [100][3];
103 float cordCarTest [100][3];
104
105 /* Generate some random Cartesian coordinates */
106 rand_m1_1((float*) cordCar, 100*3);
107
108 /* rad */
109 cart2sph((float*) cordCar, 100, SAF_FALSE, (float*) cordSph);
110 sph2cart((float*) cordSph, 100, SAF_FALSE, (float*) cordCarTest);
111 for (int i=0; i<100; i++)
112 for (int j=0; j<3; j++)
113 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, cordCar[i][j], cordCarTest[i][j]);
114
115 /* deg */
116 cart2sph((float*) cordCar, 100, SAF_TRUE, (float*) cordSph);
117 sph2cart((float*) cordSph, 100, SAF_TRUE, (float*) cordCarTest);
118 for (int i=0; i<100; i++)
119 for (int j=0; j<3; j++)
120 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, cordCar[i][j], cordCarTest[i][j]);
121}
122
123void test__delaunaynd(void){
124 int nMesh;
125 int* mesh;
126
127 /* Not really a unit test... You have to copy the mesh indices into e.g. Matlab, plot, and see... */
128
129 /* 2D 3 points */
130 float three_xy[3][2] =
131 { {7.0f,7.0f},{2.0f,7.0f},{2.0f,1.0f} };
132 mesh = NULL;
133 delaunaynd((float*)three_xy, 3, 2/*nDims*/, &mesh, &nMesh);
134 free(mesh);
135
136 /* 2D 4 points */
137 float four_xy[4][2] =
138 { {7.0f,7.0f},{2.0f,7.0f},{2.0f,1.0f},{7.0f,1.0f} };
139 mesh = NULL;
140 delaunaynd((float*)four_xy, 4, 2/*nDims*/, &mesh, &nMesh);
141 free(mesh);
142
143 /* 2D Square */
144 float square_xy[26][2] =
145 { {-1.0f,-1.0f},{-1.0f,-0.5f},{-1.0f,0.0f},{-1.0f,0.5f},{-1.0f,1.0f},{-0.5f,-1.0f},{-0.5f,-0.5f},{-0.5f,0.0f},{-0.5f,0.5f},
146 {-0.5f,1.0f},{0.0f,-1.0f},{0.0f,-0.5f},{0.0f,0.0f},{0.0f,0.5f},{0.0f,1.0f},{0.5f,-1.0f},
147 {0.5f,-0.5f},{0.5f,0.0f},{0.5f,0.5f},{0.5f,1.0f},{1.0f,-1.0f},{1.0f,-0.5f},
148 {1.0f,0.0f},{1.0f,0.5f},{1.0f,1.0f},{0.0f,0.0f} };
149 mesh = NULL;
150 delaunaynd((float*)square_xy, 26, 2/*nDims*/, &mesh, &nMesh);
151 free(mesh);
152
153 /* 3D Cube */
154 float cube_xyz[8][3] =
155 { {-1.0f,-1.0f,-1.0f},{-1.0f,1.0f,-1.0f},{1.0f,-1.0f,-1.0f},{1.0f,1.0f,-1.0f},
156 {-1.0f,-1.0f,1.0f}, {-1.0f,1.0f,1.0f}, {1.0f,-1.0f,1.0f}, {1.0f,1.0f,1.0f} };
157 mesh = NULL;
158 delaunaynd((float*)cube_xyz, 8, 3/*nDims*/, &mesh, &nMesh);
159 free(mesh);
160
161 /* 3D Cube with a point in the centre */
162 float cube_xyz2[9][3] =
163 { {-1.0f,-1.0f,-1.0f},{-1.0f,1.0f,-1.0f},{1.0f,-1.0f,-1.0f},{1.0f,1.0f,-1.0f},
164 {-1.0f,-1.0f,1.0f}, {-1.0f,1.0f,1.0f}, {1.0f,-1.0f,1.0f}, {1.0f,1.0f,1.0f}, {0.0f,0.0f,0.0f} };
165 mesh = NULL;
166 delaunaynd((float*)cube_xyz2, 9, 3/*nDims*/, &mesh, &nMesh);
167 free(mesh);
168}
169
170void test__quaternion(void){
171 int i, j;
172 float norm;
173 float rot[3][3], rot2[3][3], residual[9], ypr[3], test_ypr[3];
174 quaternion_data Q, Q1, Q2;
175
176 for(i=0; i<1000; i++){
177 /* Randomise the quaternion values */
178 rand_m1_1(Q.Q, 4);
179
180 /* Normalise to make it valid */
181 norm = L2_norm(Q.Q, 4);
182 Q.w /= norm;
183 Q.x /= norm;
184 Q.y /= norm;
185 Q.z /= norm;
186 /* Q.w = 0; Q.x = 0.0000563298236; Q.y = 0.947490811; Q.z = -0.319783032; // Problem case! */
187
188 /* Convert to rotation matrix, then back, then to rotation matrix again */
189 quaternion2rotationMatrix(&Q, rot);
190 rotationMatrix2quaternion(rot, &Q1);
191 quaternion2rotationMatrix(&Q1, rot2);
192
193 /* Ensure that the difference between them is 0 */
194 utility_svvsub((float*)rot, (float*)rot2, 9, residual);
195 for(j=0; j<9; j++)
196 TEST_ASSERT_TRUE(fabsf(residual[j])<1e-3f);
197
198 /* Testing that quaternion2euler() and euler2Quaternion() are invertable */
199 quaternion2euler(&Q1, 1, EULER_ROTATION_YAW_PITCH_ROLL, &ypr[0], &ypr[1], &ypr[2]);
200 euler2Quaternion(ypr[0], ypr[1], ypr[2], 1, EULER_ROTATION_YAW_PITCH_ROLL, &Q2);
201 quaternion2euler(&Q2, 1, EULER_ROTATION_YAW_PITCH_ROLL, &test_ypr[0], &test_ypr[1], &test_ypr[2]);
202 for(j=0; j<3; j++)
203 TEST_ASSERT_TRUE(fabsf(test_ypr[j]-ypr[j])<1e-2f);
204 }
205}
206
207void test__saf_stft_50pc_overlap(void){
208 int frame, winsize, hopsize, nFrames, ch, i, nBands, nTimesSlots, band;
209 void* hSTFT;
210 float** insig, **outsig, **inframe, **outframe;
211 float_complex*** inspec, ***outspec;
212
213 /* prep */
214 const float acceptedTolerance = 0.000001f;
215 const int fs = 48000;
216 const int signalLength = 1*fs;
217 const int framesize = 512;
218 const int nCHin = 62;
219 const int nCHout = 64;
220 insig = (float**)malloc2d(nCHin,signalLength,sizeof(float)); /* One second long */
221 outsig = (float**)malloc2d(nCHout,signalLength,sizeof(float));
222 inframe = (float**)malloc2d(nCHin,framesize,sizeof(float));
223 outframe = (float**)malloc2d(nCHout,framesize,sizeof(float));
224 rand_m1_1(FLATTEN2D(insig), nCHin*signalLength); /* populate with random numbers */
225
226 /* Set-up STFT for 50% overlapping */
227 winsize = 128;
228 hopsize = winsize/2;
229 nBands = winsize+1;
230 nTimesSlots = framesize/hopsize;
231 inspec = (float_complex***)malloc3d(nBands, nCHin, nTimesSlots, sizeof(float_complex));
232 outspec = (float_complex***)malloc3d(nBands, nCHout, nTimesSlots, sizeof(float_complex));
233 saf_stft_create(&hSTFT, winsize, hopsize, nCHin, nCHout, SAF_STFT_BANDS_CH_TIME);
234 saf_stft_channelChange(hSTFT, 123, 7); /* messing about */
235 saf_stft_flushBuffers(hSTFT); /* messing about */
236 saf_stft_channelChange(hSTFT, nCHin, nCHout); /* change back */
237
238 /* Pass insig through STFT, block-wise processing */
239 nFrames = (int)((float)signalLength/(float)framesize);
240 for(frame = 0; frame<nFrames; frame++){
241 /* Forward */
242 for(ch=0; ch<nCHin; ch++)
243 memcpy(inframe[ch], &insig[ch][frame*framesize], framesize*sizeof(float));
244 saf_stft_forward(hSTFT, inframe, framesize, inspec);
245
246 /* Copy first channel of inspec to all outspec channels */
247 for(band=0; band<nBands; band++)
248 for(ch=0; ch<nCHout; ch++)
249 memcpy(outspec[band][ch], inspec[band][0], nTimesSlots*sizeof(float_complex));
250
251 /* Backward */
252 saf_stft_backward(hSTFT, outspec, framesize, outframe);
253 for(ch=0; ch<nCHout; ch++)
254 memcpy(&outsig[ch][frame*framesize], outframe[ch], framesize*sizeof(float));
255 }
256
257 /* Check that input==output (given some numerical precision) */
258 for(i=0; i<signalLength-framesize; i++)
259 TEST_ASSERT_TRUE( fabsf(insig[0][i] - outsig[0][i+hopsize]) <= acceptedTolerance );
260
261 /* Clean-up */
262 saf_stft_destroy(&hSTFT);
263 free(insig);
264 free(outsig);
265 free(inframe);
266 free(outframe);
267 free(inspec);
268 free(outspec);
269}
270
271void test__saf_stft_LTI(void){
272 int frame, winsize, hopsize, nFrames, ch, i, nBands, nTimesSlots, band;
273 void* hSTFT;
274 float** insig, **outsig, **inframe, **outframe;
275 float_complex*** inspec, ***outspec;
276
277 /* prep */
278 const float acceptedTolerance = 0.000001f;
279 const int fs = 48000;
280 const int framesize = 128;
281 const int nCHin = 62;
282 const int nCHout = 64;
283 insig = (float**)malloc2d(nCHin,fs,sizeof(float)); /* One second long */
284 outsig = (float**)malloc2d(nCHout,fs,sizeof(float));
285 inframe = (float**)malloc2d(nCHin,framesize,sizeof(float));
286 outframe = (float**)malloc2d(nCHout,framesize,sizeof(float));
287 rand_m1_1(FLATTEN2D(insig), nCHin*fs); /* populate with random numbers */
288
289 /* Set-up STFT suitable for LTI filtering applications */
290 winsize = hopsize = 128;
291 nBands = winsize+1;
292 nTimesSlots = framesize/hopsize;
293 inspec = (float_complex***)malloc3d(nBands, nCHin, nTimesSlots, sizeof(float_complex));
294 outspec = (float_complex***)malloc3d(nBands, nCHout, nTimesSlots, sizeof(float_complex));
295 saf_stft_create(&hSTFT, winsize, hopsize, nCHin, nCHout, SAF_STFT_BANDS_CH_TIME);
296
297 /* Pass insig through STFT, block-wise processing */
298 nFrames = (int)((float)fs/(float)framesize);
299 for(frame = 0; frame<nFrames; frame++){
300 /* Forward */
301 for(ch=0; ch<nCHin; ch++)
302 memcpy(inframe[ch], &insig[ch][frame*framesize], framesize*sizeof(float));
303 saf_stft_forward(hSTFT, inframe, framesize, inspec);
304
305 /* Copy first channel of inspec to all outspec channels */
306 for(band=0; band<nBands; band++)
307 for(ch=0; ch<nCHout; ch++)
308 memcpy(outspec[band][ch], inspec[band][0], nTimesSlots*sizeof(float_complex));
309
310 /* Backward */
311 saf_stft_backward(hSTFT, outspec, framesize, outframe);
312 for(ch=0; ch<nCHout; ch++)
313 memcpy(&outsig[ch][frame*framesize], outframe[ch], framesize*sizeof(float));
314 }
315
316 /* Check that input==output (given some numerical precision) */
317 for(i=0; i<fs-framesize; i++)
318 TEST_ASSERT_TRUE( fabsf(insig[0][i] - outsig[63][i]) <= acceptedTolerance );
319
320 /* Clean-up */
321 saf_stft_destroy(&hSTFT);
322 free(insig);
323 free(outsig);
324 free(inframe);
325 free(outframe);
326 free(inspec);
327 free(outspec);
328}
329
330void test__saf_matrixConv(void){
331 int i, frame;
332 float** inputTD, **outputTD, **inputFrameTD, **outputFrameTD;
333 float*** filters;
334 void* hMatrixConv;
335
336 /* config */
337 const int signalLength = 48000;
338 const int hostBlockSize = 2048;
339 const int filterLength = 512;
340 const int nInputs = 32;
341 const int nOutputs = 40;
342
343 /* prep */
344 inputTD = (float**)malloc2d(nInputs, signalLength, sizeof(float));
345 outputTD = (float**)malloc2d(nOutputs, signalLength, sizeof(float));
346 inputFrameTD = (float**)malloc2d(nInputs, hostBlockSize, sizeof(float));
347 outputFrameTD = (float**)calloc2d(nOutputs, hostBlockSize, sizeof(float));
348 filters = (float***)malloc3d(nOutputs, nInputs, filterLength, sizeof(float));
349 rand_m1_1(FLATTEN3D(filters), nOutputs*nInputs*filterLength);
350 rand_m1_1(FLATTEN2D(inputTD), nInputs*signalLength);
351 saf_matrixConv_create(&hMatrixConv, hostBlockSize, FLATTEN3D(filters), filterLength,
352 nInputs, nOutputs, SAF_TRUE);
353
354 /* Apply */
355 for(frame = 0; frame<(int)signalLength/hostBlockSize; frame++){
356 for(i = 0; i<nInputs; i++)
357 memcpy(inputFrameTD[i], &inputTD[i][frame*hostBlockSize], hostBlockSize*sizeof(float));
358
359 saf_matrixConv_apply(hMatrixConv, FLATTEN2D(inputFrameTD), FLATTEN2D(outputFrameTD));
360
361 for(i = 0; i<nOutputs; i++)
362 memcpy(&outputTD[i][frame*hostBlockSize], outputFrameTD[i], hostBlockSize*sizeof(float));
363 }
364
365 /* Clean-up */
366 free(inputTD);
367 free(outputTD);
368 free(inputFrameTD);
369 free(outputFrameTD);
370 free(filters);
371 saf_matrixConv_destroy(&hMatrixConv);
372}
373
374void test__saf_rfft(void){
375 int i, j, N;
376 float* x_td, *test;
377 float_complex* x_fd;
378 void *hFFT;
379
380 /* Config */
381 const float acceptedTolerance = 0.00001f;
382 const int fftSizesToTest[24] =
383 {16,256,512,1024,2048,4096,8192,16384,32768,65536,1048576, /* 2^x */
384 80,160,320,640,1280,240,480,960,1920,3840,7680,15360,30720 }; /* non-2^x, (but still supported by vDSP) */
385
386 /* Loop over the different FFT sizes */
387 for (i=0; i<24; i++){
388 N = fftSizesToTest[i];
389
390 /* prep */
391 x_td = malloc1d(N*sizeof(float));
392 test = malloc1d(N*sizeof(float));
393 x_fd = malloc1d((N/2+1)*sizeof(float_complex));
394 rand_m1_1(x_td, N); /* populate with random numbers */
395 saf_rfft_create(&hFFT, N);
396
397 /* forward and backward transform */
398 saf_rfft_forward(hFFT, x_td, x_fd);
399 saf_rfft_backward(hFFT, x_fd, test);
400
401 /* Check that, x_td==test */
402 for(j=0; j<N; j++)
403 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, x_td[j], test[j]);
404
405 /* clean-up */
406 saf_rfft_destroy(&hFFT);
407 free(x_fd);
408 free(x_td);
409 free(test);
410 }
411}
412
413void test__saf_fft(void){
414 int i, j, N;
415 float_complex* x_td, *test;
416 float_complex* x_fd;
417 void *hFFT;
418
419 /* Config */
420 const float acceptedTolerance = 0.00001f;
421 const int fftSizesToTest[24] =
422 {16,256,512,1024,2048,4096,8192,16384,32768,65536,1048576, /* 2^x */
423 80,160,320,640,1280,240,480,960,1920,3840,7680,15360,30720 }; /* non-2^x, (but still supported by vDSP) */
424
425 /* Loop over the different FFT sizes */
426 for (i=0; i<24; i++){
427 N = fftSizesToTest[i];
428
429 /* prep */
430 x_td = malloc1d(N*sizeof(float_complex));
431 test = malloc1d(N*sizeof(float_complex));
432 x_fd = malloc1d(N*sizeof(float_complex));
433 rand_m1_1((float*)x_td, N*2); /* populate with random numbers */
434 saf_fft_create(&hFFT, N);
435
436 /* forward and backward transform */
437 saf_fft_forward(hFFT, x_td, x_fd);
438 saf_fft_backward(hFFT, x_fd, test);
439
440 /* Check that, x_td==test */
441 for(j=0; j<N; j++){
442 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, crealf(x_td[j]), crealf(test[j]));
443 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, cimagf(x_td[j]), cimagf(test[j]));
444 }
445
446 /* clean-up */
447 saf_fft_destroy(&hFFT);
448 free(x_fd);
449 free(x_td);
450 free(test);
451 }
452}
453
454void test__qmf(void){
455 int frame, nFrames, ch, i, nBands, procDelay, band, nHops;
456 void* hQMF;
457 float* freqVector;
458 float** insig, **outsig, **inframe, **outframe;
459 float_complex*** inspec, ***outspec;
460
461 /* prep */
462 const float acceptedTolerance = 0.01f;
463 const int fs = 48000;
464 const int signalLength = 1*fs;
465 const int framesize = 512;
466 const int hopsize = 128;
467 const int nCHin = 60;
468 const int hybridMode = 1;
469 const int nCHout = 64;
470 insig = (float**)malloc2d(nCHin,signalLength,sizeof(float)); /* One second long */
471 outsig = (float**)malloc2d(nCHout,signalLength,sizeof(float));
472 inframe = (float**)malloc2d(nCHin,framesize,sizeof(float));
473 outframe = (float**)malloc2d(nCHout,framesize,sizeof(float));
474 rand_m1_1(FLATTEN2D(insig), nCHin*signalLength); /* populate with random numbers */
475
476 /* Set-up */
477 nHops = framesize/hopsize;
478 qmf_create(&hQMF, nCHin, nCHout, hopsize, hybridMode, QMF_BANDS_CH_TIME);
479 procDelay = qmf_getProcDelay(hQMF);
480 nBands = qmf_getNBands(hQMF);
481 freqVector = malloc1d(nBands*sizeof(float));
482 qmf_getCentreFreqs(hQMF, (float)fs, nBands, freqVector);
483 inspec = (float_complex***)malloc3d(nBands, nCHin, nHops, sizeof(float_complex));
484 outspec = (float_complex***)malloc3d(nBands, nCHout, nHops, sizeof(float_complex));
485
486 /* Pass insig through the QMF filterbank, block-wise processing */
487 nFrames = (int)((float)signalLength/(float)framesize);
488 for(frame = 0; frame<nFrames; frame++){
489 /* QMF Analysis */
490 for(ch=0; ch<nCHin; ch++)
491 memcpy(inframe[ch], &insig[ch][frame*framesize], framesize*sizeof(float));
492 qmf_analysis(hQMF, inframe, framesize, inspec);
493
494 /* Copy first channel of inspec to all outspec channels */
495 for(band=0; band<nBands; band++)
496 for(ch=0; ch<nCHout; ch++)
497 memcpy(outspec[band][ch], inspec[band][0], nHops*sizeof(float_complex));
498
499 /* QMF Synthesis */
500 qmf_synthesis(hQMF, outspec, framesize, outframe);
501 for(ch=0; ch<nCHout; ch++)
502 memcpy(&outsig[ch][frame*framesize], outframe[ch], framesize*sizeof(float));
503 }
504
505 /* Check that input==output (given some numerical precision) - channel 0 */
506 for(i=0; i<signalLength-procDelay-framesize; i++)
507 TEST_ASSERT_TRUE( fabsf(insig[0][i] - outsig[0][i+procDelay]) <= acceptedTolerance );
508
509 /* Clean-up */
510 qmf_destroy(&hQMF);
511 free(insig);
512 free(outsig);
513 free(inframe);
514 free(outframe);
515 free(inspec);
516 free(outspec);
517 free(freqVector);
518}
519
520void test__smb_pitchShifter(void){
521 float* inputData, *outputData;
522 void* hPS, *hFFT;
523 float frequency;
524 int i, smbLatency, ind;
525
526 /* Config */
527 const int sampleRate = 48000;
528 const int FFTsize = 8192;
529 const int osfactor = 4;
530 const int nSamples = 8*FFTsize;
531
532 /* prep */
533 smb_pitchShift_create(&hPS, 1, FFTsize, osfactor, (float)sampleRate);
534 inputData = malloc1d(nSamples*sizeof(float));
535 outputData = calloc1d(nSamples,sizeof(float));
536 frequency = (float)sampleRate/8.0f;
537 for(i=0; i<nSamples; i++) /* sine tone at quarter Nyquist: */
538 inputData[i] = sinf(2.0f * SAF_PI * (float)i * frequency/(float)sampleRate);
539 smbLatency = FFTsize - (FFTsize/osfactor);
540
541 /* Pitch shift down one octave */
542 smb_pitchShift_apply(hPS, 0.5, nSamples, inputData, outputData);
543
544 /* Take FFT, the bin with the highest energy should correspond to 1/8 Nyquist */
545 float_complex* out_fft; // [nSamples / 2 + 1];
546 out_fft = malloc1d((nSamples / 2 + 1) * sizeof(float_complex));
547 saf_rfft_create(&hFFT, nSamples);
548 saf_rfft_forward(hFFT, outputData, out_fft);
549 utility_cimaxv(out_fft, nSamples/2+1, &ind);
550 TEST_ASSERT_TRUE(ind == nSamples/16);
551
552 /* clean-up */
553 smb_pitchShift_destroy(&hPS);
554 saf_rfft_destroy(&hFFT);
555 free(inputData);
556 free(outputData);
557 free(out_fft);
558}
559
560void test__sortf(void){
561 float* values;
562 int* sortedIdx;
563 int i;
564
565 /* Config */
566 const int numValues = 10000;
567
568 /* Prep */
569 sortedIdx = malloc1d(numValues*sizeof(int));
570 values = malloc1d(numValues*sizeof(float));
571 rand_m1_1(values, numValues); /* populate with random numbers */
572
573 /* Sort in accending order */
574 sortf(values, NULL, sortedIdx, numValues, 0);
575
576 /* Check that the next value is either the same or greater than current value */
577 for(i=0; i<numValues-1; i++)
578 TEST_ASSERT_TRUE(values[sortedIdx[i]]<=values[sortedIdx[i+1]]);
579
580 /* Sort in decending order */
581 sortf(values, NULL, sortedIdx, numValues, 1);
582
583 /* Check that the next value is either the same or less than current value */
584 for(i=0; i<numValues-1; i++)
585 TEST_ASSERT_TRUE(values[sortedIdx[i]]>=values[sortedIdx[i+1]]);
586
587 /* clean-up */
588 free(values);
589 free(sortedIdx);
590}
591
592void test__sortz(void){
593 int i;
594 const int N = 36;
595 double_complex vals[36] ={
596 cmplx(1.0, 1.0), cmplx(7.0, 1.0), cmplx(10.0, 5.0),
597 cmplx(12.0, 4.0), cmplx(4.0, 4.0), cmplx(8.0, 0.0),
598 cmplx(10.0, -1.0), cmplx(7.0, 5.0), cmplx(7.0, 2.0),
599 cmplx(5.0, 1.0), cmplx(4.0, -1.0), cmplx(32.0, 3.0),
600 cmplx(32.0, 32.5), cmplx(25.0, 0.0), cmplx(2.0, -2.0),
601 cmplx(7.0, -2.0), cmplx(1.0, -1.0), cmplx(12.0, -1.0),
602 cmplx(2.0, -1.0), cmplx(4.0, 2.0), cmplx(10.0, 6.0),
603 cmplx(5.0, 2.0), cmplx(32.0, 1.5), cmplx(7.0, -10.0),
604 cmplx(1.0, -1.5), cmplx(4.0, 25.0), cmplx(3.0, 2.0),
605 cmplx(1.0, 4.5), cmplx(10.0, 5.0), cmplx(10.0, 2.0),
606 cmplx(10.0, -3.5), cmplx(30.0, -10.0), cmplx(7.0, -12.0),
607 cmplx(1.0, -13.5), cmplx(12.0, -12.0), cmplx(32.0, 23.0)
608 };
609 double_complex sorted_vals[36];
610
611 /* Sort assending order */
612 sortz(vals, sorted_vals, N, 0);
613
614 /* Check that the next real(value) is either the same or greater than current real(value) */
615 for(i=0; i<N-1; i++)
616 TEST_ASSERT_TRUE(creal(sorted_vals[i])<=creal(sorted_vals[i+1]));
617
618 /* Check that if the next real(value) is the same as the current real(value), then
619 * the next imag(value) is greater that the current image(value)*/
620 for(i=0; i<N-1; i++)
621 if(fabs(creal(sorted_vals[i])-creal(sorted_vals[i+1])) < 0.00001 )
622 TEST_ASSERT_TRUE(cimag(sorted_vals[i])<=cimag(sorted_vals[i+1]));
623
624 /* Sort decending order */
625 sortz(vals, sorted_vals, N, 1);
626
627 /* Check that the next real(value) is either the same or smaller than current real(value) */
628 for(i=0; i<N-1; i++)
629 TEST_ASSERT_TRUE(creal(sorted_vals[i])>=creal(sorted_vals[i+1]));
630
631 /* Check that if the next real(value) is the same as the current real(value), then
632 * the next imag(value) is smaller that the current image(value)*/
633 for(i=0; i<N-1; i++)
634 if(fabs(creal(sorted_vals[i])-creal(sorted_vals[i+1])) < 0.00001 )
635 TEST_ASSERT_TRUE(cimag(sorted_vals[i])>=cimag(sorted_vals[i+1]));
636}
637
638void test__cmplxPairUp(void){
639 int i;
640 const int N = 36;
641 double_complex vals[36] ={
642 cmplx(1.0, 1.0), cmplx(7.0, 1.0), cmplx(10.0, 5.0),
643 cmplx(12.0, 4.0), cmplx(4.0, 4.0), cmplx(8.0, 0.0),
644 cmplx(10.0, -1.0), cmplx(7.0, 5.0), cmplx(7.0, 2.0),
645 cmplx(5.0, 1.0), cmplx(4.0, -1.0), cmplx(32.0, 3.0),
646 cmplx(32.0, 32.5), cmplx(25.0, 0.0), cmplx(2.0, -2.0),
647 cmplx(7.0, -2.0), cmplx(1.0, -1.0), cmplx(12.0, -1.0),
648 cmplx(2.0, -1.0), cmplx(4.0, 2.0), cmplx(10.0, 6.0),
649 cmplx(5.0, 0.0), cmplx(32.0, 1.5), cmplx(7.0, -10.0),
650 cmplx(1.0, -1.5), cmplx(4.0, 25.0), cmplx(3.0, 2.0),
651 cmplx(1.0, 0.0), cmplx(10.0, 5.0), cmplx(10.0, 2.0),
652 cmplx(10.0, -3.5), cmplx(30.0, -10.0), cmplx(7.0, -12.0),
653 cmplx(1.0, -13.5), cmplx(12.0, -12.0), cmplx(32.0, 23.0)
654 };
655 double_complex sorted_vals[36];
656
657 /* Sort assending order */
658 cmplxPairUp(vals, sorted_vals, N);
659
660 /* Check that the next real(value) is either the same or greater than current real(value),
661 * Ignoring purely real numbers */
662 for(i=0; i<N-1; i++)
663 if( !(fabs(cimag(sorted_vals[i])) < 0.0001) && !(fabs(cimag(sorted_vals[i+1])) < 0.0001) )
664 TEST_ASSERT_TRUE(creal(sorted_vals[i])<=creal(sorted_vals[i+1]));
665
666 /* Check that the next real(value) is either the same or greater than current real(value),
667 * Only considering purely real numbers */
668 for(i=0; i<N-1; i++)
669 if( (fabs(cimag(sorted_vals[i])) < 0.0001) && (fabs(cimag(sorted_vals[i+1])) < 0.0001) )
670 TEST_ASSERT_TRUE(creal(sorted_vals[i])<=creal(sorted_vals[i+1]));
671
672 /* Check that if the next real(value) is the same as the current real(value), then
673 * the next imag(value) is greater that the current image(value)
674 * Ignoring purely real numbers */
675 for(i=0; i<N-1; i++)
676 if(fabs(creal(sorted_vals[i])-creal(sorted_vals[i+1])) < 0.00001 )
677 if( !(fabs(cimag(sorted_vals[i])) < 0.0001) && !(fabs(cimag(sorted_vals[i+1])) < 0.0001) )
678 TEST_ASSERT_TRUE(cimag(sorted_vals[i])<=cimag(sorted_vals[i+1]));
679}
680
681void test__getVoronoiWeights(void){
682 int i, it, td, nDirs;
683 float* dirs_deg, *weights;
684 float sum, tmp, scale;
685
686 /* Config */
687 const float acceptedTolerance = 0.01f;
688 const int nIterations = 100;
689
690 /* Loop over T-designs */
691 for(td=2; td<21; td++){
692 dirs_deg = (float*)__HANDLES_Tdesign_dirs_deg[td];
693 nDirs = __Tdesign_nPoints_per_degree[td];
694
695 /* Compute weights */
696 weights = malloc1d(nDirs*sizeof(float));
697 getVoronoiWeights(dirs_deg, nDirs, 0, weights);
698
699 /* Assert that they sum to 4PI */
700 sum = sumf(weights, nDirs);
701 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 4.0f*SAF_PI, sum);
702
703 /* Due to the uniform arrangement, all the weights should be the same */
704 for(i=1; i<nDirs; i++)
705 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, weights[0], weights[i]);
706
707 /* clean-up */
708 free(weights);
709 }
710
711 /* Loop over some random arrangement of points */
712 for(it=0; it<nIterations; it++){
713 rand_0_1(&tmp, 1);
714 nDirs = (int)(tmp*190.0f + 10.0f); /* random number between 10..200 */
715
716 /* Random dirs (-180..180 azi, -180..180 elev) */
717 dirs_deg = malloc1d(nDirs*2*sizeof(float));
718 rand_m1_1(dirs_deg, nDirs*2);
719 scale = 180.0f;
720 utility_svsmul(dirs_deg, &scale, nDirs*2, dirs_deg);
721
722 /* Compute weights */
723 weights = malloc1d(nDirs*sizeof(float));
724 getVoronoiWeights(dirs_deg, nDirs, 0, weights);
725
726 /* Assert that they sum to 4PI */
727 sum = sumf(weights, nDirs);
728 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, 4.0f*SAF_PI, sum);
729
730 /* clean-up */
731 free(dirs_deg);
732 free(weights);
733 }
734}
735
736void test__unique_i(void){
737 int i, nUnique;
738 int* uniqueVals;
739 int* uniqueInds;
740
741 /* test1 */
742 int input[6] = {1, 2, 2, 10, 11, 12};
743 int uniqueVals_ref[5] = {1, 2, 10, 11, 12};
744 int uniqueInds_ref[5] = {0, 2, 3, 4, 5};
745 unique_i(input, 6, &uniqueVals, &uniqueInds, &nUnique);
746 TEST_ASSERT_EQUAL(5, nUnique);
747 for(i=0; i<nUnique; i++){
748 TEST_ASSERT_EQUAL(uniqueVals_ref[i], uniqueVals[i]);
749 TEST_ASSERT_EQUAL(uniqueInds_ref[i], uniqueInds[i]);
750 }
751 free(uniqueVals);
752 free(uniqueInds);
753
754 /* test2 */
755 int input2[12] = {1, 10, 1, 3, 1, 3, 4, 7, 8, 10, 10, 2};
756 int uniqueVals_ref2[7] = {1, 3, 4, 7, 8, 10, 2};
757 int uniqueInds_ref2[7] = {4, 5, 6, 7, 8, 10, 11};
758 unique_i(input2, 12, &uniqueVals, &uniqueInds, &nUnique);
759 TEST_ASSERT_EQUAL(7, nUnique);
760 for(i=0; i<nUnique; i++){
761 TEST_ASSERT_EQUAL(uniqueVals_ref2[i], uniqueVals[i]);
762 TEST_ASSERT_EQUAL(uniqueInds_ref2[i], uniqueInds[i]);
763 }
764 free(uniqueVals);
765 free(uniqueInds);
766}
767
768void test__latticeDecorrelator(void){
769 int c, band, nBands, idx, hopIdx, i;
770 void* hDecor, *hSTFT;
771 float icc, tmp, tmp2;
772 float* freqVector;
773 float** inputTimeDomainData, **outputTimeDomainData, **tempHop;
774 float_complex*** inTFframe, ***outTFframe;
775
776 /* config */
777 const float acceptedICC = 0.05f;
778 const int nCH = 24;
779 const int nTestHops = 800;
780 const int hopSize = 128;
781 const int procDelay = hopSize*12 + 12; // from afSTFT + decor
782 const int lSig = nTestHops*hopSize+procDelay;
783 const float fs = 48e3f;
784 nBands = hopSize+5;
785
786 /* audio buffers */
787 inputTimeDomainData = (float**) calloc2d(1, lSig, sizeof(float));
788 outputTimeDomainData = (float**) calloc2d(nCH, lSig, sizeof(float));
789 inTFframe = (float_complex***)malloc3d(nBands, nCH, 1, sizeof(float_complex));
790 outTFframe = (float_complex***)malloc3d(nBands, nCH, 1, sizeof(float_complex));
791 tempHop = (float**) malloc2d(nCH, hopSize, sizeof(float));
792
793 /* Initialise afSTFT and input data */
794 afSTFT_create(&hSTFT, 1, nCH, hopSize, 0, 1, AFSTFT_BANDS_CH_TIME);
795 rand_m1_1(FLATTEN2D(inputTimeDomainData), 1*lSig); /* populate with random numbers */
796 freqVector = malloc1d(nBands*sizeof(float));
797 afSTFT_getCentreFreqs(hSTFT, fs, nBands, freqVector);
798
799 /* setup decorrelator */
800 int orders[4] = {20, 15, 6, 6}; /* 20th order up to 700Hz, 15th->2.4kHz, 6th->4kHz, 3rd->12kHz, NONE(only delays)->Nyquist */
801 //float freqCutoffs[4] = {600.0f, 2.6e3f, 4.5e3f, 12e3f};
802 float freqCutoffs[4] = {900.0f, 6.8e3f, 12e3f, 24e3f};
803 const int maxDelay = 12;
804 latticeDecorrelator_create(&hDecor, fs, hopSize, freqVector, nBands, nCH, orders, freqCutoffs, 4, maxDelay, 0, 0.75f);
805
806 /* Processing loop */
807 idx = 0;
808 hopIdx = 0;
809 while(idx<lSig-hopSize*2){
810 for(c=0; c<1; c++)
811 memcpy(tempHop[c], &(inputTimeDomainData[c][hopIdx*hopSize]), hopSize*sizeof(float));
812
813 /* forward TF transform, and replicate to all channels */
814 afSTFT_forward(hSTFT, tempHop, hopSize, inTFframe);
815 for(band=0; band<nBands; band++)
816 for(i=1; i<nCH;i++)
817 inTFframe[band][i][0] = inTFframe[band][0][0];
818
819 /* decorrelate */
820 latticeDecorrelator_apply(hDecor, inTFframe, 1, outTFframe);
821
822 /* backward TF transform */
823 afSTFT_backward(hSTFT, outTFframe, hopSize, tempHop);
824
825 /* Copy frame to output TD buffer */
826 for(c=0; c<nCH; c++)
827 memcpy(&(outputTimeDomainData[c][hopIdx*hopSize]), tempHop[c], hopSize*sizeof(float));
828 idx+=hopSize;
829 hopIdx++;
830 }
831
832 /* Compensate for processing delay, and check that the inter-channel correlation
833 * coefficient is below the accepted threshold (ideally 0, if fully
834 * decorrelated...) */
835 for(c=0; c<nCH; c++){
836 utility_svvdot(inputTimeDomainData[0], &outputTimeDomainData[c][procDelay], (lSig-procDelay), &icc);
837 utility_svvdot(inputTimeDomainData[0], inputTimeDomainData[0], (lSig-procDelay), &tmp);
838 utility_svvdot(&outputTimeDomainData[c][procDelay], &outputTimeDomainData[c][procDelay], (lSig-procDelay), &tmp2);
839
840 icc = icc/sqrtf(tmp*tmp2); /* normalise */
841 TEST_ASSERT_TRUE(fabsf(icc)<acceptedICC);
842 }
843#if 0
844 /* Check for mutually decorrelated channels... */
845 int c2;
846 for(c=0; c<nCH; c++){
847 for(c2=0; c2<nCH; c2++){
848 utility_svvdot(&outputTimeDomainData[c][procDelay], &outputTimeDomainData[c2][procDelay], (lSig-procDelay), &icc);
849 utility_svvdot(&outputTimeDomainData[c2][procDelay], &outputTimeDomainData[c2][procDelay], (lSig-procDelay), &tmp);
850 utility_svvdot(&outputTimeDomainData[c][procDelay], &outputTimeDomainData[c][procDelay], (lSig-procDelay), &tmp2);
851
852 if (c!=c2){
853 icc = icc/sqrtf(tmp*tmp2); /* normalise */
854 TEST_ASSERT_TRUE(fabsf(icc)<acceptedICC);
855 }
856 }
857 }
858#endif
859
860 /* Clean-up */
861 latticeDecorrelator_destroy(&hDecor);
862 free(inTFframe);
863 free(outTFframe);
864 free(tempHop);
865 afSTFT_destroy(&hSTFT);
866 free(inputTimeDomainData);
867 free(outputTimeDomainData);
868}
869
870void test__butterCoeffs(void){
871 int i;
872 float fs, cutoff_freq, cutoff_freq2;
873 int order;
874
875 /* Config */
876 const double acceptedTolerance = 0.00001f;
877
878 /* 1st order Low-pass filter */
879 fs = 48e3f;
880 cutoff_freq = 3000.0f;
881 order = 1;
882 double a_test1[2], b_test1[2];
883 butterCoeffs(BUTTER_FILTER_LPF, order, cutoff_freq, 0.0f, fs, (double*)b_test1, (double*)a_test1);
884 const double a_ref1[2] = {1,-0.668178637919299};
885 const double b_ref1[2] = {0.165910681040351,0.165910681040351};
886 for(i=0; i<2; i++){ /* Compare with the values given by Matlab's butter function */
887 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a_test1[i], a_ref1[i]);
888 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b_test1[i], b_ref1[i]);
889 }
890
891 /* 2nd order Low-pass filter */
892 fs = 48e3f;
893 cutoff_freq = 12000.0f;
894 order = 2;
895 double a_test2[3], b_test2[3];
896 butterCoeffs(BUTTER_FILTER_LPF, order, cutoff_freq, 0.0f, fs, (double*)b_test2, (double*)a_test2);
897 const double a_ref2[3] = {1.0,-2.22044604925031e-16,0.171572875253810};
898 const double b_ref2[3] = {0.292893218813452,0.585786437626905,0.292893218813452};
899 for(i=0; i<3; i++){ /* Compare with the values given by Matlab's butter function */
900 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a_test2[i], a_ref2[i]);
901 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b_test2[i], b_ref2[i]);
902 }
903
904 /* 3rd order Low-pass filter */
905 fs = 48e3f;
906 cutoff_freq = 200.0f;
907 order = 3;
908 double a_test3[4], b_test3[4];
909 butterCoeffs(BUTTER_FILTER_LPF, order, cutoff_freq, 0.0f, fs, (double*)b_test3, (double*)a_test3);
910 const double a_ref3[4] = {1.0,-2.94764161678340,2.89664496645376,-0.948985866903327};
911 const double b_ref3[4] = {2.18534587909103e-06,6.55603763727308e-06,6.55603763727308e-06,2.18534587909103e-06};
912 for(i=0; i<4; i++){ /* Compare with the values given by Matlab's butter function */
913 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a_test3[i], a_ref3[i]);
914 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b_test3[i], b_ref3[i]);
915 }
916
917 /* 6th order Low-pass filter */
918 fs = 48e3f;
919 cutoff_freq = 1e3f;
920 order = 6;
921 double a_test4[7], b_test4[7];
922 butterCoeffs(BUTTER_FILTER_LPF, order, cutoff_freq, 0.0f, fs, (double*)b_test4, (double*)a_test4);
923 const double a_ref4[7] = {1,-5.49431292177096,12.5978414666894,-15.4285267903275,10.6436770055305,-3.92144696766748,0.602772146971300};
924 const double b_ref4[7] = {6.15535184628202e-08,3.69321110776921e-07,9.23302776942303e-07,1.23107036925640e-06,9.23302776942303e-07,3.69321110776921e-07,6.15535184628202e-08};
925 for(i=0; i<7; i++){ /* Compare with the values given by Matlab's butter function */
926 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a_test4[i], a_ref4[i]);
927 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b_test4[i], b_ref4[i]);
928 }
929
930 /* 3rd order High-pass filter */
931 fs = 48e3f;
932 cutoff_freq = 3000.0f;
933 order = 3;
934 double a_test5[4], b_test5[4];
935 butterCoeffs(BUTTER_FILTER_HPF, order, cutoff_freq, 0.0f, fs, (double*)b_test5, (double*)a_test5);
936 const double a_ref5[4] = {1,-2.21916861831167,1.71511783003340,-0.453545933365530};
937 const double b_ref5[4] = {0.673479047713825,-2.02043714314147,2.02043714314147,-0.673479047713825};
938 for(i=0; i<4; i++){ /* Compare with the values given by Matlab's butter function */
939 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a_test5[i], a_ref5[i]);
940 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b_test5[i], b_ref5[i]);
941 }
942
943 /* 4th order High-pass filter */
944 fs = 48e3f;
945 cutoff_freq = 100.0;
946 order = 4;
947 double a_test6[5], b_test6[5];
948 butterCoeffs(BUTTER_FILTER_HPF, order, cutoff_freq, 0.0f, fs, (double*)b_test6, (double*)a_test6);
949 const double a_ref6[5] = {1.0,-3.96579438007005,5.89796693861409,-3.89854491737242,0.966372387692057};
950 const double b_ref6[5] = {0.983042413984288,-3.93216965593715,5.89825448390573,-3.93216965593715,0.983042413984288};
951 for(i=0; i<5; i++){ /* Compare with the values given by Matlab's butter function */
952 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a_test6[i], a_ref6[i]);
953 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b_test6[i], b_ref6[i]);
954 }
955
956 /* 2nd order Band-pass filter */
957 fs = 48e3f;
958 cutoff_freq = 100.0;
959 cutoff_freq2 = 400.0;
960 order = 2;
961 double a_test7[5], b_test7[5];
962 butterCoeffs(BUTTER_FILTER_BPF, order, cutoff_freq, cutoff_freq2, fs, (double*)b_test7, (double*)a_test7);
963 const double a_ref7[5] = {1.0,-3.94312581006024,5.83226704209421,-3.83511871130750,0.945977936232284};
964 const double b_ref7[5] = {0.000375069616051004,0.0,-0.000750139232102008,0.0,0.000375069616051004};
965 for(i=0; i<5; i++){ /* Compare with the values given by Matlab's butter function */
966 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a_test7[i], a_ref7[i]);
967 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b_test7[i], b_ref7[i]);
968 }
969
970 /* 3rd order Band-stop filter */
971 fs = 48e3f;
972 cutoff_freq = 240.0;
973 cutoff_freq2 = 1600.0;
974 order = 3;
975 double a_test9[7], b_test9[7];
976 butterCoeffs(BUTTER_FILTER_BSF, order, cutoff_freq, cutoff_freq2, fs, (double*)b_test9, (double*)a_test9);
977 const double a_ref9[7] = {1,-5.62580309774365,13.2124846784594,-16.5822627287366,11.7304049556188,-4.43493124452282,0.700107676775329};
978 const double b_ref9[7] = {0.836724592951539,-5.00379660039217,12.4847741945760,-16.6354041344203,12.4847741945760,-5.00379660039217,0.836724592951539};
979 for(i=0; i<7; i++){ /* Compare with the values given by Matlab's butter function */
980 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a_test9[i], a_ref9[i]);
981 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b_test9[i], b_ref9[i]);
982 }
983}
984
985/* Test evalIIRTransferFunction()
986 * Coefficients for the first 7 tests below are taken from the butterworth tests
987 * above, so can be compared against the mag/phase results given by MATLAB's
988 * freqz function. The 8th test loop evaluates results for 1st order shelving
989 * filters with coefficients generated by the DVF filter functions. The last
990 * test loop runs those same tests, but with the floating point version of
991 * evalIIRTransferFunctionf() (valid for low order filters)
992 */
993void test__evalIIRTransferFunction(void) {
994 int i, nCoeffs;
995 float phase_ref, mag_ref, mag_ref_db;
996
997 /* Config */
998 const double phaseTolerance = 0.0174533f * 5.f; /* ~ 1 degree * mul */
999 const double magToleranceDb = 0.1f; /* tolerance in dB, for a target magnitude of 0dB */
1000 const float errScale = 2.f / 120.f; /* tolerance grows for lower dB target: toleranceLevel/atLevel.
1001 * e.g. 2/120 = 2dB tolerance for -120 target dB */
1002 const int nFreqs = 10;
1003 const float freqs[10] = {147.21423f,270.49564f,411.40091f,687.90202f,1395.3072f,2024.3936f,3696.9416f,6784.4745f,9798.67f,17594.058f};
1004 float mag[10], phase[10];
1005 float fs = 44.1e3f;
1006
1007 /* evalIIRTransferFunction(): coeffs of type double */
1008
1009 /* Test 1 * 1st order Low-pass filter */
1010 nCoeffs = 2;
1011 const double a_t1[2] = {1,-0.6681786};
1012 const double b_t1[2] = {0.1659107,0.1659107};
1013 const double mag_ref1[10] = {0.99861294,0.99533929,0.98931332,0.97092312,0.89393904,0.80765661,0.59366549,0.35440473,0.23070415,0.06521195};
1014 const double phase_ref1[10] = {-0.052676048,-0.096585083,-0.14632679,-0.24173907,-0.46473796,-0.63062928,-0.93519006,-1.2085189,-1.337995,-1.5055381};
1015 evalIIRTransferFunction((double*)b_t1, (double*)a_t1, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase);
1016 for(i = 0; i < nFreqs; i++) {
1017 phase_ref = (float)phase_ref1[i];
1018 mag_ref = (float)mag_ref1[i];
1019 mag_ref_db = 20.f*log10f(mag_ref);
1020 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1021 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1022 }
1023
1024 /* Test 2 * 2nd order Low-pass filter */
1025 nCoeffs = 3;
1026 const double a_t2[3] = {1,-0,0.1715729};
1027 const double b_t2[3] = {0.2928932,0.5857864,0.2928932};
1028 const float mag_ref2[10] = {0.99999991f,0.99999985f,0.99999955f,0.99999702f,0.99995046f,0.99977761f,0.99736787f,0.96409579f,0.81776268f,0.1073164f};
1029 const float phase_ref2[10] = {-0.014832279f,-0.027258003f,-0.041470589f,-0.069414188f,-0.14150066f,-0.20679977f,-0.39012494f,-0.7974393f,-1.3261562f,-2.6614056f};
1030 evalIIRTransferFunction((double*)b_t2, (double*)a_t2, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase);
1031 for(i = 0; i < nFreqs; i++) {
1032 phase_ref = phase_ref2[i];
1033 mag_ref = mag_ref2[i];
1034 mag_ref_db = 20.f*log10f(mag_ref);
1035 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1036 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1037 }
1038
1039 /* Test 3 * 3rd order Low-pass filter */
1040 nCoeffs = 4;
1041 const double a_t3[4] = {1,-2.9476416,2.896645,-0.9489859};
1042 const double b_t3[4] = {2.2e-06,6.6e-06,6.6e-06,2.2e-06};
1043 const float mag_ref3[10] = {0.8954383f,0.3011618f,0.0892913f,0.0191409f,0.0022769f,0.0007374f,0.0001152f,1.56e-05f,3.8e-06f,1e-07f};
1044 const float phase_ref3[10] = {-1.8249735f,3.0678618f,2.4995092f,2.1114704f,1.8340934f,1.751328f,1.6679375f,1.6206872f,1.6020054f,1.579398f};
1045 evalIIRTransferFunction((double*)b_t3, (double*)a_t3, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase);
1046 for(i = 0; i < nFreqs; i++) {
1047 phase_ref = phase_ref3[i];
1048 mag_ref = mag_ref3[i];
1049 mag_ref_db = 20.f*log10f(mag_ref);
1050 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1051 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1052 }
1053
1054 /* Test 4 * 6th order Low-pass filter */
1055 nCoeffs = 7;
1056 const double a_t4[7] = {1,-5.49431292177096,12.5978414666894,-15.4285267903275,10.6436770055305,-3.92144696766748,0.6027721469713};
1057 const double b_t4[7] = {6.15535184628202e-08,3.69321110776921e-07,9.23302776942303e-07,1.2310703692564e-06,9.23302776942303e-07,3.69321110776921e-07,6.15535184628202e-08};
1058 const float mag_ref4[10] = {0.9999999834907868f,0.9999997831836054f,0.9999679556869572f,0.9849426248859378f,0.08033081621985783f,0.008452216914022819f,0.0002063542729228268f,3.793812554381118e-06f,2.274031694371124e-07f,9.970589432354785e-11f};
1059 const float phase_ref4[10] = {-0.6201852189230334f,-1.148525513374147f,-1.774695369143539f,3.109543344373707f,-0.4296773811384472f,-1.349824316530828f,-2.195405632723407f,-2.65814688739603f,-2.839508904295157f,-3.058387834019209f};
1060 evalIIRTransferFunction((double*)b_t4, (double*)a_t4, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase);
1061 for(i = 0; i < nFreqs; i++) {
1062 float phase_ref = phase_ref4[i];
1063 float mag_ref = mag_ref4[i];
1064 float mag_ref_db = 20.f*log10f(mag_ref);
1065 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1066 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1067 }
1068
1069 /* Test 5 * 3rd order High-pass filter */
1070 nCoeffs = 4;
1071 const double a_t5[4] = {1,-2.2191686,1.7151178,-0.4535459};
1072 const double b_t5[4] = {0.673479,-2.0204371,2.0204371,-0.673479};
1073 const float mag_ref5[10] = {0.0001466f,0.0009096f,0.0032014f,0.0149875f,0.125037f,0.362653f,0.927991f,0.9985214f,0.9999112f,0.9999999f};
1074 const float phase_ref5[10] = {-1.6762949f,-1.7648759f,-1.866651f,-2.0692621f,-2.6256366f,3.0800183f,1.6530258f,0.7789431f,0.4789307f,0.1307956f};
1075 evalIIRTransferFunction((double*)b_t5, (double*)a_t5, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase);
1076 for(i = 0; i < nFreqs; i++) {
1077 float phase_ref = phase_ref5[i];
1078 float mag_ref = mag_ref5[i];
1079 float mag_ref_db = 20.f*log10f(mag_ref);
1080 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1081 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1082 }
1083
1084 /* Test 6 * 4th order High-pass filter
1085 * 400 Hz cut (differs from butterworth test above) */
1086 nCoeffs = 5;
1087 const double a_t6[5] = {1,-3.863184622426,5.598835456747838,-3.607752453919942,0.872108645089876};
1088 const double b_t6[5] = {4.3909323578772e-07,1.75637294315089e-06,2.63455941472633e-06,1.75637294315089e-06,4.3909323578772e-07};
1089 const float mag_ref6[10] = {0.9996691528983467f,0.9595570109649983f,0.5370184819357747f,0.08100263003740536f,0.004753436194609436f,0.001057169058757887f,8.896712774518116e-05f,6.197328265811134e-06f,9.491865964914827e-07f,5.478157027512644e-09f};
1090 const float phase_ref6[10] = {-1.072517623166929f,-2.13344694428915f,2.732267641095127f,1.462991201859678f,0.6929733816699927f,0.4733493046075806f,0.2541184532330854f,0.130425028023503f,0.08157492611996242f,0.02248140228360206f};
1091 evalIIRTransferFunction((double*)b_t6, (double*)a_t6, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase);
1092 for(i = 0; i < nFreqs; i++) {
1093 float phase_ref = phase_ref6[i];
1094 float mag_ref = mag_ref6[i];
1095 float mag_ref_db = 20.f*log10f(mag_ref);
1096 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1097 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1098 }
1099
1100 /* Test 7 * 2nd order Band-pass filter */
1101 nCoeffs = 5;
1102 const double a_t7[5] = {1,-3.9431258,5.832267,-3.8351187,0.9459779};
1103 const double b_t7[5] = {0.0003751,0,-0.0007501,0,0.0003751};
1104 const float mag_ref7[10] = {0.7829909f,0.9051549f,0.5636772f,0.1816557f,0.0400635f,0.0185759f,0.0053305f,0.0014022f,0.0005484f,4.16e-05f};
1105 const float phase_ref7[10] = {0.4017825f,-0.7852502f,-1.8127451f,-2.4983166f,-2.8544848f,-2.9475768f,-3.0381483f,-3.0886103f,-3.1084696f,-3.1324667f};
1106 evalIIRTransferFunction((double*)b_t7, (double*)a_t7, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase); for(i=0; i<nFreqs; i++){
1107 float phase_ref = phase_ref7[i];
1108 float mag_ref = mag_ref7[i];
1109 float mag_ref_db = 20.f*log10f(mag_ref);
1110 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1111 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1112 }
1113
1114 /* Test the response for 12 settings used by the DVF filters
1115 * (1st order sheving filter) */
1116 const int nTest = 12;
1117 nCoeffs = 2;
1118 double as_dvf[12][2] = { /* nTest x nCoeffs */
1119 {1,-0.95864619},{1,-0.96599375},{1,-0.9648155},{1,-0.84969467},{1,-0.93327999},{1,-0.95974372},
1120 {1,-0.83460338},{1,-0.74744027},{1,-0.67445272},{1,-0.76911048},{1,-0.64266857},{1,-0.54043336}
1121 };
1122 double bs_dvf[12][2] = { /* nTest x nCoeffs */
1123 {8.0841171,-7.5217374},{2.1404179,-2.0439339},{1.3371466,-1.2837154},{0.73353449,-0.52570938},
1124 {1.1312827,-1.0324739},{1.1334695,-1.0826754},{0.18784397,-0.098191093},{0.43493823,-0.24012268},
1125 {0.72850398,-0.42469436},{0.1577158,-0.075691788},{0.34545179,-0.16259809},{0.60618525,-0.27852873}
1126 };
1127 const float mags_dvf[12][10] = { /* nTest x nFreqs */
1128 {12.68472f,11.390262f,10.245287f,9.081507f,8.2881923f,8.1239126f,8.0139907f,7.9799796f,7.9724969f,7.9680408f},
1129 {2.6652868f,2.469902f,2.3321513f,2.2179604f,2.1523794f,2.139901f,2.1317514f,2.1292619f,2.1287161f,2.1283915f},
1130 {1.473629f,1.4224625f,1.3865239f,1.35693f,1.3400517f,1.3368519f,1.3347643f,1.3341269f,1.3339872f,1.3339041f},
1131 {1.3740782f,1.3545086f,1.3209933f,1.2349497f,1.0204791f,0.89934391f,0.76424296f,0.70524874f,0.69060103f,0.68154193f},
1132 {1.4538524f,1.4034324f,1.3412611f,1.2506542f,1.1632856f,1.1414561f,1.125968f,1.1210229f,1.1199249f,1.1192693f},
1133 {1.2358328f,1.2022487f,1.1755623f,1.151321f,1.1364575f,1.1335478f,1.1316331f,1.1310458f,1.130917f,1.1308403f},
1134 {0.53871826f,0.53107297f,0.51772931f,0.48194542f,0.3814555f,0.3150888f,0.22673701f,0.17897735f,0.16548424f,0.15666687f},
1135 {0.76984932f,0.76629998f,0.75986062f,0.74092839f,0.6717326f,0.60914565f,0.49873472f,0.42504312f,0.40260097f,0.38760994f},
1136 {0.93261062f,0.93115748f,0.92849149f,0.92042692f,0.88786149f,0.85374491f,0.78101307f,0.72261208f,0.7032242f,0.68987066f},
1137 {0.3542684f,0.3519694f,0.34781956f,0.33577027f,0.29342253f,0.25693916f,0.1950882f,0.15408756f,0.14134236f,0.13268952f},
1138 {0.51134341f,0.51045389f,0.50881513f,0.50380352f,0.48269293f,0.45893406f,0.4014785f,0.34644275f,0.32577048f,0.31064565f},
1139 {0.71281398f,0.71244818f,0.7117705f,0.70966741f,0.70027544f,0.68857561f,0.65428371f,0.61109598f,0.59151358f,0.57580457f}
1140 };
1141 const float phases_dvf[12][10] = { /* nTest x nFreqs */
1142 {-0.17782001f,-0.24874011f,-0.2637155f,-0.22716735f,-0.13811864f,-0.098858451f,-0.054719219f,-0.028348151f,-0.01776687f,-0.0049023852f},
1143 {-0.11818797f,-0.14315719f,-0.13342746f,-0.10042039f,-0.055483624f,-0.038914731f,-0.021247104f,-0.010960723f,-0.006863078f,-0.001892663f},
1144 {-0.054685172f,-0.064793725f,-0.059294582f,-0.043867626f,-0.023978567f,-0.016782288f,-0.0091501707f,-0.0047182622f,-0.0029540703f,-0.00081461202f},
1145 {-0.064892969f,-0.11661773f,-0.17043929f,-0.25417001f,-0.34353751f,-0.33903683f,-0.25655254f,-0.15106232f,-0.097685198f,-0.027478557f},
1146 {-0.069273584f,-0.10993957f,-0.13346817f,-0.13654806f,-0.096248577f,-0.071345378f,-0.040468966f,-0.0211288f,-0.013264997f,-0.0036639447f},
1147 {-0.042920762f,-0.054369797f,-0.052367865f,-0.040534678f,-0.022767208f,-0.01601877f,-0.0087641883f,-0.0045240297f,-0.0028331222f,-0.00078136769f},
1148 {-0.082361693f,-0.14918816f,-0.22112994f,-0.34301241f,-0.52772513f,-0.5813414f,-0.53771111f,-0.3682489f,-0.25049777f,-0.073003569f},
1149 {-0.036111073f,-0.06587829f,-0.098882791f,-0.15880356f,-0.2712657f,-0.32156729f,-0.32729224f,-0.23402028f,-0.16071005f,-0.047082247f},
1150 {-0.014103031f,-0.025780252f,-0.03883755f,-0.063048578f,-0.11207771f,-0.13775677f,-0.14908283f,-0.11045335f,-0.07655027f,-0.022550556f},
1151 {-0.050350716f,-0.09181969f,-0.13772611f,-0.22079247f,-0.37596653f,-0.44682619f,-0.46554397f,-0.34629266f,-0.24235026f,-0.072094835f},
1152 {-0.019043671f,-0.03486822f,-0.052685779f,-0.086317621f,-0.15955383f,-0.2051296f,-0.24897791f,-0.20818672f,-0.15155625f,-0.046353162f},
1153 {-0.006828925f,-0.012518705f,-0.018958244f,-0.031275577f,-0.059563883f,-0.079317458f,-0.10560439f,-0.097624085f,-0.0740598f,-0.023377524f}
1154 };
1155 for(int t = 0; t < nTest; t++) {
1156 float phase_ref, mag_ref, mag_ref_db;
1157 double* pAs = &as_dvf[t][0];
1158 double* pBs = &bs_dvf[t][0];
1159 evalIIRTransferFunction(pBs, pAs, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase);
1160 for(i = 0; i < nFreqs; i++){
1161 phase_ref = phases_dvf[t][i];
1162 mag_ref = mags_dvf[t][i];
1163 mag_ref_db = 20.f*log10f(mag_ref);
1164 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1165 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1166 }
1167 }
1168 /* using the same parameters as above, test evalIIRTransferFunctionf():
1169 * coeffs of type float */
1170 float as_dvf_f[12][2] = { /* nTest x nCoeffs */
1171 {1.f,-0.95864619f},{1.f,-0.96599375f},{1.f,-0.9648155f},{1.f,-0.84969467f},{1.f,-0.93327999f},{1.f,-0.95974372f},
1172 {1.f,-0.83460338f},{1.f,-0.74744027f},{1.f,-0.67445272f},{1.f,-0.76911048f},{1.f,-0.64266857f},{1.f,-0.54043336f}
1173 };
1174 float bs_dvf_f[12][2] = { /* nTest x nCoeffs */
1175 {8.0841171f,-7.5217374f},{2.1404179f,-2.0439339f},{1.3371466f,-1.2837154f},{0.73353449f,-0.52570938f},
1176 {1.1312827f,-1.0324739f},{1.1334695f,-1.0826754f},{0.18784397f,-0.098191093f},{0.43493823f,-0.24012268f},
1177 {0.72850398f,-0.42469436f},{0.1577158f,-0.075691788f},{0.34545179f,-0.16259809f},{0.60618525f,-0.27852873f}
1178 };
1179 for(int t = 0; t < nTest; t++) {
1180 float* pAs = &as_dvf_f[t][0];
1181 float* pBs = &bs_dvf_f[t][0];
1182 evalIIRTransferFunctionf(pBs, pAs, nCoeffs, (float*)freqs, nFreqs, fs, 0, (float*)mag, (float*)phase);
1183 for(i = 0; i < nFreqs; i++){
1184 phase_ref = phases_dvf[t][i];
1185 mag_ref = mags_dvf[t][i];
1186 mag_ref_db = 20.f*log10f(mag_ref);
1187 TEST_ASSERT_FLOAT_WITHIN(magToleranceDb + errScale*fabsf(mag_ref_db), mag_ref_db, 20*log10(mag[i]));
1188 TEST_ASSERT_FLOAT_WITHIN(phaseTolerance, phase_ref, phase[i]);
1189 }
1190 }
1191
1192}
1193
1194
1195void test__faf_IIRFilterbank(void){
1196 void* hFaF;
1197 int i, band;
1198 float* inSig, *outSig;
1199 float** outFrame, **outSig_bands;
1200 void* hFFT;
1201 float_complex* insig_fft, *outsig_fft;
1202
1203 /* Config */
1204 const float acceptedTolerance_dB = 0.5f;
1205 const int signalLength = 256;
1206 const int frameSize = 256;//16;
1207 float fs = 48e3;
1208 int order = 3;
1209 float fc[6] = {176.776695296637f, 353.553390593274f, 707.106781186547f, 1414.21356237309f, 2828.42712474619f, 5656.85424949238f};
1210 inSig = malloc1d(signalLength * sizeof(float));
1211 outSig_bands = (float**)malloc2d(7, signalLength, sizeof(float));
1212 outSig = calloc1d(signalLength, sizeof(float));
1213
1214 insig_fft = malloc1d((signalLength / 2 + 1) * sizeof(float_complex));
1215 outsig_fft = malloc1d((signalLength / 2 + 1) * sizeof(float_complex));
1216
1217 /* Impulse */
1218 memset(inSig, 0, signalLength*sizeof(float));
1219 inSig[0] = 1.0f;
1220
1221 /* Pass impulse through filterbank */
1222 outFrame = (float**)malloc2d(7, frameSize, sizeof(float));
1223 faf_IIRFilterbank_create(&hFaF, order, (float*)fc, 6, fs, frameSize);
1224 for(i=0; i< signalLength/frameSize; i++){
1225 faf_IIRFilterbank_apply(hFaF, &inSig[i*frameSize], outFrame, frameSize);
1226 for(band=0; band<7; band++)
1227 memcpy(&outSig_bands[band][i*frameSize], outFrame[band], frameSize*sizeof(float));
1228 }
1229 faf_IIRFilterbank_destroy(&hFaF);
1230
1231 /* Sum the individual bands */
1232 for(band=0; band<7; band++)
1233 cblas_saxpy(signalLength, 1.0f, outSig_bands[band], 1, outSig, 1);
1234
1235 /* Check that the magnitude difference between input and output is below 0.5dB */
1236 saf_rfft_create(&hFFT, signalLength);
1237 saf_rfft_forward(hFFT, inSig, insig_fft);
1238 saf_rfft_forward(hFFT, outSig, outsig_fft);
1239 for(i=0; i<signalLength/2+1; i++)
1240 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance_dB, 0.0f, 20.0f * log10f(cabsf( ccdivf(outsig_fft[i],insig_fft[i]) )));
1241
1242 /* Now the same thing, but for 1st order */
1243 order = 1;
1244 faf_IIRFilterbank_create(&hFaF, order, (float*)fc, 6, fs, frameSize);
1245 for(i=0; i< signalLength/frameSize; i++){
1246 faf_IIRFilterbank_apply(hFaF, &inSig[i*frameSize], outFrame, frameSize);
1247 for(band=0; band<7; band++)
1248 memcpy(&outSig_bands[band][i*frameSize], outFrame[band], frameSize*sizeof(float));
1249 }
1250 faf_IIRFilterbank_destroy(&hFaF);
1251 memset(outSig, 0, signalLength*sizeof(float));
1252 for(band=0; band<7; band++)
1253 cblas_saxpy(signalLength, 1.0f, outSig_bands[band], 1, outSig, 1);
1254 saf_rfft_forward(hFFT, outSig, outsig_fft);
1255 for(i=0; i<signalLength/2+1; i++)
1256 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance_dB, 0.0f, 20.0f * log10f(cabsf(ccdivf(outsig_fft[i], insig_fft[i]))));
1257
1258 /* clean-up */
1259 saf_rfft_destroy(&hFFT);
1260 free(outFrame);
1261 free(inSig);
1262 free(outSig_bands);
1263 free(outSig);
1264 free(insig_fft);
1265 free(outsig_fft);
1266}
1267
1268void test__gexpm(void){
1269 int i, j;
1270 float outM[6][6];
1271
1272 /* prep */
1273 const float acceptedTolerance = 0.0001f;
1274 const float inM[6][6] = {
1275 {-0.376858200853762f,0.656790634216694f,0.124479178614046f,-0.334752428307223f,1.50745241578235f,0.0290651989052969f},
1276 {0.608382058262806f,0.581930485432986f,3.23135406998058f,-0.712003744668929f,-1.33872571354702f,-0.334742482743222f},
1277 {-0.795741418256672f,0.690709474622409f,0.620971281129248f,1.38749471231620f,0.897245329198841f,-0.0693670166113321f},
1278 {0.179789913109994f,-1.06135084902804f,-1.10032635271188f,0.612441344250358f,-2.43213807790664f,-0.479265889956047f},
1279 {-0.277441781278754f,-0.0732116130293688f,-0.572551795688137f,1.02024767389969f,0.167385894565923f,1.45210312619277f},
1280 {-0.205305770089918f,-1.59783032780633f,1.08539265129120f,0.460057585947626f,-1.02420974042838f,1.04117461500218f}
1281 };
1282 const float outM_ref[6][6] = {
1283 {0.385163650730121f,0.0865151585709784f,0.898406722231524f,0.877640791713973f,0.435244824708340f,0.888866982998854f},
1284 {-0.664938511314777f,5.02943129352875f,8.24444951891833f,2.23840978101979f,-0.942669833528886f,-2.38535530623266f},
1285 {-0.388189314743059f,0.429308537172675f,1.13870842882926f,1.60875776611798f,-1.44249911796405f,-1.51822150286392f},
1286 {1.05630187656688f,0.256606570814868f,-2.42701873560847f,-1.42372526577009f,-0.335273289873574f,-1.94362909671742f},
1287 {0.0261470437116839f,-3.03329326250434f,-3.50207776203591f,0.412043775125377f,-0.536000387729306f,1.61801775548557f},
1288 {-0.292024827617294f,-4.31537192033477f,-3.99160103133879f,0.312499067924889f,-1.46924802440347f,1.98522802303672f}
1289 };
1290
1291 /* Compute matrix exponential */
1292 gexpm((float*)inM, 6, 0, (float*)outM);
1293
1294 /* Check that output of SAF's gexpm, is similar to Matlab's expm: */
1295 for(i=0; i<6; i++)
1296 for(j=0; j<6; j++)
1297 TEST_ASSERT_TRUE( fabsf(outM[i][j] - outM_ref[i][j]) <= acceptedTolerance );
1298}
1299
1300/* Target values are generated by MATLAB functions in `generate_coeffs_for_plugin_tests.m`
1301 * (corresponding function names are noted above each data set), which is not included in this
1302 * repository but are in the author's development repository `nearfield_rangeextrapolation`.
1303 */
1304void test__dvf_calcDVFShelfParams(void){
1305
1306 /* setup */
1307 const float acceptedTolerance = 0.00001f;
1308 const float acceptedTolerance_fc = 0.1f;
1309 const int nTheta = 19;
1310 const int nRho = 5;
1311 const float rho[5] = {1.150000f,1.250000f,1.570000f,2.381000f,3.990000f};
1312 const float g0_ref[5][19] = { /* testRhoCoeffs_g_0 */
1313 {22.670282f,17.717752f,11.902597f,7.906282f,4.720884f,2.217061f,0.134088f,-1.613730f,-3.095960f,-5.279052f,-5.433653f,-6.342905f,-7.107677f,-7.744796f,-8.236084f,-8.613662f,-8.864276f,-9.065870f,-9.089385f},
1314 {18.295933f,15.510441f,11.452312f,7.951083f,4.997235f,2.609075f,0.592613f,-1.107964f,-2.557504f,-4.547256f,-4.853912f,-5.750024f,-6.504702f,-7.133244f,-7.621092f,-7.993574f,-8.244015f,-8.438287f,-8.467470f},
1315 {11.937032f,11.093339f,9.245757f,7.118216f,4.990070f,3.083402f,1.371444f,-0.121838f,-1.427296f,-2.979742f,-3.542803f,-4.381065f,-5.091220f,-5.683427f,-6.149122f,-6.508598f,-6.748356f,-6.923465f,-6.961620f},
1316 {6.676990f,6.451424f,5.818377f,4.924700f,3.861979f,2.760683f,1.662668f,0.629080f,-0.327831f,-1.328149f,-1.970549f,-2.649238f,-3.234743f,-3.727775f,-4.122829f,-4.433178f,-4.640902f,-4.783351f,-4.823625f},
1317 {3.628860f,3.534311f,3.298166f,2.922799f,2.438587f,1.888286f,1.296135f,0.698518f,0.112899f,-0.473626f,-0.960644f,-1.428032f,-1.841763f,-2.196404f,-2.487131f,-2.717121f,-2.873915f,-2.978235f,-3.010937f}
1318 };
1319 const float gInf_ref[5][19] = { /* testRhoCoeffs_g_inf */
1320 {-4.643651f,-4.225287f,-4.134752f,-4.386332f,-5.244711f,-6.439307f,-7.659091f,-8.887172f,-10.004796f,-10.694171f,-11.190476f,-10.876569f,-10.140292f,-9.913242f,-9.411469f,-8.981807f,-8.723677f,-8.529900f,-8.574359f},
1321 {-4.128221f,-3.834507f,-3.575000f,-3.637788f,-4.278932f,-5.310000f,-6.609705f,-7.815000f,-8.925450f,-9.646588f,-10.000000f,-9.784733f,-9.301643f,-8.862963f,-8.370815f,-7.953778f,-7.693305f,-7.500645f,-7.518260f},
1322 {-3.094135f,-2.963709f,-2.721834f,-2.573043f,-2.793627f,-3.414652f,-4.403297f,-5.518539f,-6.578461f,-7.332562f,-7.702192f,-7.582977f,-7.376856f,-6.745349f,-6.279895f,-5.891862f,-5.636418f,-5.456323f,-5.437006f},
1323 {-1.937289f,-1.889079f,-1.765709f,-1.620800f,-1.598110f,-1.815613f,-2.314443f,-3.041183f,-3.857777f,-4.533446f,-4.931544f,-4.962571f,-4.717069f,-4.357935f,-3.971281f,-3.646312f,-3.422461f,-3.269044f,-3.231471f},
1324 {-1.126412f,-1.103440f,-1.049199f,-0.969714f,-0.917898f,-0.962176f,-1.182409f,-1.566237f,-2.065834f,-2.552771f,-2.884909f,-2.977707f,-2.811758f,-2.629199f,-2.355800f,-2.118920f,-1.949860f,-1.834291f,-1.800638f}
1325 };
1326 const float fc_ref[5][19] = { /* testRhoCoeffs_f_c */
1327 {525.636204f,409.078426f,427.552571f,936.671283f,1635.128987f,2622.394035f,3167.199181f,3899.649293f,4176.703569f,4361.226917f,4634.448076f,4516.401848f,4567.834168f,4685.234222f,4908.786495f,4966.258562f,4936.982049f,4927.963688f,5210.861482f},
1328 {410.072475f,389.319631f,398.844102f,613.394238f,1116.223303f,2095.651724f,2847.557763f,3726.141143f,4080.406901f,4304.960791f,4463.911798f,4449.375495f,4501.166349f,4623.582375f,4757.884246f,4911.093999f,4935.074404f,4940.266143f,5155.085794f},
1329 {358.247441f,352.931439f,352.752741f,402.566754f,590.733021f,1127.131294f,2007.589994f,3160.896502f,3808.131027f,4155.246718f,4336.853155f,4375.553567f,4406.656373f,4543.636509f,4649.878140f,4849.374583f,4974.986343f,5006.214905f,5164.504029f},
1330 {318.842699f,318.199637f,315.776327f,326.423309f,364.498469f,500.548368f,980.626776f,2174.301881f,3296.777215f,3904.656864f,4203.152454f,4329.347194f,4338.652755f,4492.976051f,4579.879128f,4849.678327f,5052.801340f,5116.753611f,5267.402018f},
1331 {297.454930f,297.570719f,296.701047f,300.362959f,308.255747f,342.596563f,509.934390f,1379.970914f,2702.827191f,3646.599635f,4078.866661f,4301.570222f,4303.807248f,4472.223890f,4535.654099f,4855.399825f,5119.558700f,5210.329993f,5380.750972f}
1332 };
1333 /* params to test */
1334 float g0, gInf, fc;
1335
1336 for(int ri = 0; ri < nRho; ri++){
1337 float r = rho[ri];
1338 for(int ti = 0; ti < nTheta; ti++){
1339 calcDVFShelfParams(ti, r, &g0, &gInf, &fc);
1340 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, g0_ref[ri][ti], g0);
1341 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, gInf_ref[ri][ti], gInf);
1342 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance_fc, fc_ref[ri][ti], fc);
1343 }
1344 }
1345}
1346
1347/* Parameter interpolation is implicitly checked in test__dvf_dvfShelfCoeffs. */
1348void test__dvf_interpDVFShelfParams(void){
1349
1350 /* interpDVFShelfParams() calls calcDVFShelfParams() twice to generate the
1351 * high shelf parameters for the nearests angles in the lookup table. Those
1352 * parameters are subsequently interpolated. So the success of
1353 * interpHighShelfCoeffs() relies first on the calcDVFShelfParams(),
1354 * so that should be tested first. */
1355
1356 /* setup */
1357 const float acceptedTolerance = 0.0001f;
1358 const float acceptedTolerance_frq = 0.01f;
1359 const int nTheta = 6;
1360 const float theta[6] = {0.000000f,2.300000f,47.614000f,98.600000f,166.200000f,180.000000f};
1361 const int nRho = 5;
1362 const float rho[5] = {1.150000f,1.250000f,1.570000f,2.381000f,3.990000f};
1363 const float iG0_ref[5][6] = { /* testShelfParamsInterp_iG_0 */
1364 {22.670282f,21.531200f,2.814473f,-5.412009f,-8.989264f,-9.089385f},
1365 {18.295933f,17.655270f,3.178890f,-4.810981f,-8.364464f,-8.467470f},
1366 {11.937032f,11.742982f,3.538333f,-3.463974f,-6.856924f,-6.961620f},
1367 {6.676990f,6.625110f,3.023452f,-1.880613f,-4.729220f,-4.823625f},
1368 {3.628860f,3.607114f,2.019588f,-0.892461f,-2.938593f,-3.010937f}
1369 };
1370 const float iGInf_ref[5][6] = { /* testShelfParamsInterp_iG_inf */
1371 {-4.643651f,-4.547427f,-6.154277f,-11.120993f,-8.603536f,-8.574359f},
1372 {-4.128221f,-4.060667f,-5.063987f,-9.950522f,-7.573856f,-7.518260f},
1373 {-3.094135f,-3.064137f,-3.266476f,-7.650444f,-5.524759f,-5.437006f},
1374 {-1.937289f,-1.926201f,-1.763717f,-4.875811f,-3.327342f,-3.231471f},
1375 {-1.126412f,-1.121129f,-0.951611f,-2.838410f,-1.878207f,-1.800638f}
1376 };
1377 const float iFc_ref[5][6] = { /* testShelfParamsInterp_if_c */
1378 {525.636204f,498.827915f,2386.832594f,4596.197114f,4931.390665f,5210.861482f},
1379 {410.072475f,405.299321f,1861.960103f,4441.658657f,4938.293282f,5155.085794f},
1380 {358.247441f,357.024760f,999.146666f,4311.428254f,4994.348051f,5164.504029f},
1381 {318.842699f,318.694795f,468.086862f,4161.363071f,5092.451748f,5267.402018f},
1382 {297.454930f,297.481562f,334.402844f,4018.349277f,5175.836901f,5380.750972f}
1383 };
1384 /* High shelf parameters to be tested */
1385 float iG0, iGInf, iFc;
1386
1387 for(int ri = 0; ri < nRho; ri++){
1388 float r = rho[ri];
1389 for(int ti = 0; ti < nTheta; ti++){
1390 interpDVFShelfParams(theta[ti], r, &iG0, &iGInf, &iFc);
1391 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, iG0, iG0_ref[ri][ti]);
1392 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, iGInf, iGInf_ref[ri][ti]);
1393 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance_frq, iFc, iFc_ref[ri][ti]);
1394 }
1395 }
1396}
1397
1398void test__dvf_dvfShelfCoeffs(void){
1399 /* setup */
1400 const float acceptedTolerance = 0.00001f;
1401 const int nTheta = 6;
1402 const float theta[6] = {0.000000f,2.300000f,47.614000f,98.600000f,166.200000f,180.000000f};
1403 const int nRho = 5;
1404 const float rho[5] = {1.150000f,1.250000f,1.570000f,2.381000f,3.990000f};
1405 const float fs = 44100;
1406 const float b0_ref[5][6] = { /* testIIRCoefs_b0 */
1407 {8.084117f,7.162779f,0.733534f,0.181211f,0.157716f,0.157753f},
1408 {5.162983f,4.832104f,0.847478f,0.218379f,0.188114f,0.188100f},
1409 {2.787888f,2.735502f,1.052975f,0.322169f,0.274301f,0.274359f},
1410 {1.733188f,1.725027f,1.162720f,0.508950f,0.432220f,0.432405f},
1411 {1.337147f,1.334601f,1.133469f,0.693396f,0.606185f,0.606313f}
1412 };
1413 const float b1_ref[5][6] = { /* testIIRCoefs_b1 */
1414 {-7.521737f,-6.689086f,-0.525709f,-0.092147f,-0.075692f,-0.072026f},
1415 {-4.880667f,-4.570874f,-0.654751f,-0.113974f,-0.090171f,-0.086752f},
1416 {-2.654257f,-2.604818f,-0.917824f,-0.171818f,-0.130191f,-0.126320f},
1417 {-1.659057f,-1.651278f,-1.090421f,-0.278191f,-0.201595f,-0.195404f},
1418 {-1.283715f,-1.281267f,-1.082675f,-0.387883f,-0.278529f,-0.268341f}
1419 };
1420 const float a1_ref[5][6] = { /* testIIRCoefs_a1 */
1421 {-0.958646f,-0.960287f,-0.849695f,-0.833925f,-0.769110f,-0.755889f},
1422 {-0.965649f,-0.965782f,-0.866341f,-0.818335f,-0.743436f,-0.731349f},
1423 {-0.966189f,-0.966188f,-0.910070f,-0.775971f,-0.682649f,-0.670043f},
1424 {-0.965632f,-0.965605f,-0.948954f,-0.713458f,-0.602472f,-0.587018f},
1425 {-0.964816f,-0.964791f,-0.959744f,-0.661426f,-0.540433f,-0.522001f}
1426 };
1427
1428 float b0, b1, a1; /* IIR coeffs to be tested */
1429 float iG0, iGInf, iFc;
1430
1431 for(int ri = 0; ri < nRho; ri++){
1432 float r = rho[ri];
1433 for(int ti = 0; ti < nTheta; ti++){
1434 interpDVFShelfParams(theta[ti], r, &iG0, &iGInf, &iFc);
1435 dvfShelfCoeffs(iG0, iGInf, iFc, fs, &b0, &b1, &a1);
1436 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b0, b0_ref[ri][ti]);
1437 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, b1, b1_ref[ri][ti]);
1438 TEST_ASSERT_FLOAT_WITHIN(acceptedTolerance, a1, a1_ref[ri][ti]);
1439 }
1440 }
1441}
afSTFT_create
void afSTFT_create(void **const phSTFT, int nCHin, int nCHout, int hopsize, int lowDelayMode, int hybridmode, AFSTFT_FDDATA_FORMAT format)
Creates an instance of afSTFT.
Definition afSTFTlib.c:143
afSTFT_forward
void afSTFT_forward(void *const hSTFT, float **dataTD, int framesize, float_complex ***dataFD)
Performs forward afSTFT transform.
Definition afSTFTlib.c:230
afSTFT_getCentreFreqs
void afSTFT_getCentreFreqs(void *const hSTFT, float fs, int nBands, float *freqVector)
Returns current frequency vector.
Definition afSTFTlib.c:546
afSTFT_destroy
void afSTFT_destroy(void **const phSTFT)
Destroys an instance of afSTFT.
Definition afSTFTlib.c:199
afSTFT_backward
void afSTFT_backward(void *const hSTFT, float_complex ***dataFD, int framesize, float **dataTD)
Performs backward afSTFT transform.
Definition afSTFTlib.c:348
AFSTFT_BANDS_CH_TIME
@ AFSTFT_BANDS_CH_TIME
nBands x nChannels x nTimeHops
Definition afSTFTlib.h:80
faf_IIRFilterbank_destroy
void faf_IIRFilterbank_destroy(void **phFaF)
Destroys an instance of the Favrot & Faller filterbank.
Definition saf_utility_filters.c:1241
smb_pitchShift_create
void smb_pitchShift_create(void **hSmb, int nCH, int fftFrameSize, int osamp, float sampleRate)
Creates an instance of SMB PitchShifter.
Definition saf_utility_pitch.c:124
saf_matrixConv_apply
void saf_matrixConv_apply(void *const hMC, float *inputSig, float *outputSig)
Performs the matrix convolution.
Definition saf_utility_matrixConv.c:178
saf_stft_backward
void saf_stft_backward(void *const hSTFT, float_complex ***dataFD, int framesize, float **dataTD)
Performs the backward-STFT operation for the current frame.
Definition saf_utility_fft.c:439
SAF_PI
#define SAF_PI
pi constant (single precision)
Definition saf_utilities.h:70
saf_stft_forward
void saf_stft_forward(void *const hSTFT, float **dataTD, int framesize, float_complex ***dataFD)
Performs the forward-STFT operation for the current frame.
Definition saf_utility_fft.c:367
__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
sortf
void sortf(float *in_vec, float *out_vec, int *new_idices, int len, int descendFLAG)
Sort a vector of floating-point values into ascending/decending order (optionally returning the new i...
Definition saf_utility_sort.c:153
SAF_FALSE
#define SAF_FALSE
Boolean false.
Definition saf_utilities.h:67
quaternion2rotationMatrix
void quaternion2rotationMatrix(quaternion_data *Q, float R[3][3])
Constructs a 3x3 rotation matrix based on a quaternion.
Definition saf_utility_geometry.c:90
interpDVFShelfParams
void interpDVFShelfParams(float theta, float rho, float *iG0, float *iGInf, float *iFc)
Calculate the shelving filter parameters for the Distance Variation Function filter from the source (...
Definition saf_utility_dvf.c:108
bessel_Jn
void bessel_Jn(int N, double *z, int nZ, double *J_n, double *dJ_n)
Computes the values of the (cylindrical) Bessel function of the first kind (Jn) and it's derivative (...
Definition saf_utility_bessel.c:381
qmf_getCentreFreqs
void qmf_getCentreFreqs(void *const hQMF, float fs, int nBands, float *centreFreq)
Computes the QMF/hybrid-QMF centre frequencies.
Definition saf_utility_qmf.c:628
saf_fft_backward
void saf_fft_backward(void *const hFFT, float_complex *inputFD, float_complex *outputTD)
Performs the backward-FFT operation; use for complex to complex transformations.
Definition saf_utility_fft.c:935
euler2Quaternion
void euler2Quaternion(float alpha, float beta, float gamma, int degreesFlag, EULER_ROTATION_CONVENTIONS convention, quaternion_data *Q)
Converts Euler angles to a quaternion.
Definition saf_utility_geometry.c:124
unique_i
void unique_i(int *input, int nInputs, int **uniqueVals, int **uniqueInds, int *nUnique)
Finds the unique values (and their indices) of the input vector.
Definition saf_utility_misc.c:404
getVoronoiWeights
void getVoronoiWeights(float *dirs_deg, int nDirs, int diagFLAG, float *weights)
Computes the integration weights, based on the areas of each face of the corresponding Voronoi diagra...
Definition saf_utility_geometry.c:944
utility_svvsub
void utility_svvsub(const float *a, const float *b, const int len, float *c)
Single-precision, vector-vector subtraction, i.e.
Definition saf_utility_veclib.c:842
cmplxPairUp
void cmplxPairUp(double_complex *in_vec, double_complex *out_vec, int len)
Pairs up complex numbers and sorts them in ascending order based on their real parts first,...
Definition saf_utility_sort.c:341
saf_fft_create
void saf_fft_create(void **const phFFT, int N)
Creates an instance of saf_fft; complex<->complex FFT.
Definition saf_utility_fft.c:762
qmf_getProcDelay
int qmf_getProcDelay(void *const hQMF)
Returns the processing delay in samples.
Definition saf_utility_qmf.c:610
saf_stft_channelChange
void saf_stft_channelChange(void *const hSTFT, int new_nCHin, int new_nCHout)
Changes the number of input/output channels.
Definition saf_utility_fft.c:491
smb_pitchShift_apply
void smb_pitchShift_apply(void *hSmb, float pitchShift, int frameSize, float *indata, float *outdata)
Performs pitch shifting of the input signals, while retaining the same time duration as the original ...
Definition saf_utility_pitch.c:213
qmf_synthesis
void qmf_synthesis(void *const hQMF, float_complex ***dataFD, int framesize, float **dataTD)
Performs QMF synthesis of the input frequency-domain signals.
Definition saf_utility_qmf.c:438
saf_matrixConv_create
void saf_matrixConv_create(void **const phMC, int hopSize, float *H, int length_h, int nCHin, int nCHout, int usePartFLAG)
Creates an instance of matrixConv.
Definition saf_utility_matrixConv.c:50
quaternion2euler
void quaternion2euler(quaternion_data *Q, int degreesFlag, EULER_ROTATION_CONVENTIONS convention, float *alpha, float *beta, float *gamma)
Converts a quaternion to Euler angles.
Definition saf_utility_geometry.c:164
faf_IIRFilterbank_create
void faf_IIRFilterbank_create(void **phFaF, int order, float *fc, int nCutoffFreq, float sampleRate, int maxNumSamples)
Computes a bank of IIR filter coefficients required to divide a signal into frequency bands,...
Definition saf_utility_filters.c:1047
rand_0_1
void rand_0_1(float *vector, int length)
Generates random numbers between 0 and 1 and stores them in the input vector.
Definition saf_utility_misc.c:250
utility_cimaxv
void utility_cimaxv(const float_complex *a, const int len, int *index)
Single-precision, complex, index of maximum absolute value in a vector, i.e.
Definition saf_utility_veclib.c:401
saf_rfft_create
void saf_rfft_create(void **const phFFT, int N)
Creates an instance of saf_rfft; real<->half-complex (conjugate-symmetric) FFT.
Definition saf_utility_fft.c:532
SAF_TRUE
#define SAF_TRUE
Boolean true.
Definition saf_utilities.h:64
sph2cart
void sph2cart(float *sph, int nDirs, int anglesInDegreesFLAG, float *cart)
Converts spherical coordinates to Cartesian coordinates.
Definition saf_utility_geometry.c:272
smb_pitchShift_destroy
void smb_pitchShift_destroy(void **const hSmb)
Destroys an instance of SMB PitchShifter.
Definition saf_utility_pitch.c:170
qmf_analysis
void qmf_analysis(void *const hQMF, float **dataTD, int framesize, float_complex ***dataFD)
Performs QMF analysis of the input time-domain signals.
Definition saf_utility_qmf.c:315
saf_fft_forward
void saf_fft_forward(void *const hFFT, float_complex *inputTD, float_complex *outputFD)
Performs the forward-FFT operation; use for complex to complex transformations.
Definition saf_utility_fft.c:898
gexpm
void gexpm(float *D, int sizeD, int m1, float *Y)
Numerically solves first-order, linear, homogeneous differential equation systems,...
Definition saf_utility_misc.c:512
evalIIRTransferFunctionf
void evalIIRTransferFunctionf(float *b_coeff, float *a_coeff, int nCoeffs, float *freqs, int nFreqs, float fs, int mag2dB, float *magnitude, float *phase_rad)
Computes magnitude and phase response of an IIR filter from its coefficients (floats) at user-specifi...
Definition saf_utility_filters.c:610
rotationMatrix2quaternion
void rotationMatrix2quaternion(float R[3][3], quaternion_data *Q)
Calculates the quaternion corresponding to a 3x3 rotation matrix.
Definition saf_utility_geometry.c:108
qmf_destroy
void qmf_destroy(void **const phQMF)
Destroys an instance of the qmf filterbank.
Definition saf_utility_qmf.c:272
latticeDecorrelator_apply
void latticeDecorrelator_apply(void *hDecor, float_complex ***inFrame, int nTimeSlots, float_complex ***decorFrame)
Applies the lattice all-pass-filter-based multi-channel signal decorrelator.
Definition saf_utility_decor.c:341
saf_matrixConv_destroy
void saf_matrixConv_destroy(void **const phMC)
Destroys an instance of matrixConv.
Definition saf_utility_matrixConv.c:133
saf_rfft_forward
void saf_rfft_forward(void *const hFFT, float *inputTD, float_complex *outputFD)
Performs the forward-FFT operation; use for real to complex (conjugate symmetric) transformations.
Definition saf_utility_fft.c:670
saf_rfft_backward
void saf_rfft_backward(void *const hFFT, float_complex *inputFD, float *outputTD)
Performs the backward-FFT operation; use for complex (conjugate symmetric) to real transformations.
Definition saf_utility_fft.c:715
sumf
float sumf(float *values, int nValues)
Returns the sum of all values.
Definition saf_utility_misc.c:372
saf_rfft_destroy
void saf_rfft_destroy(void **const phFFT)
Destroys an instance of saf_rfft.
Definition saf_utility_fft.c:617
qmf_getNBands
int qmf_getNBands(void *const hQMF)
Returns the number of frequency bands.
Definition saf_utility_qmf.c:619
latticeDecorrelator_create
void latticeDecorrelator_create(void **phDecor, float fs, int hopsize, float *freqVector, int nBands, int nCH, int *orders, float *freqCutoffs, int nCutoffs, int maxDelay, int lookupOffset, float enComp_coeff)
Creates an instance of the lattice all-pass-filter-based multi-channel signal decorrelator.
Definition saf_utility_decor.c:195
saf_stft_destroy
void saf_stft_destroy(void **const phSTFT)
Destroys an instance of saf_stft.
Definition saf_utility_fft.c:347
butterCoeffs
void butterCoeffs(BUTTER_FILTER_TYPES filterType, int order, float cutoff1, float cutoff2, float sampleRate, double *b_coeffs, double *a_coeffs)
Computes Butterworth IIR filter coefficients [1].
Definition saf_utility_filters.c:814
sortz
void sortz(double_complex *in_vec, double_complex *out_vec, int len, int descendFLAG)
Sort a vector of complex double floating-point values into ascending/ decending order.
Definition saf_utility_sort.c:277
dvfShelfCoeffs
void dvfShelfCoeffs(float g0, float gInf, float fc, float fs, float *b0, float *b1, float *a1)
Calculate the DVF filter coefficients from shelving filter parameters.
Definition saf_utility_dvf.c:149
__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
bessel_yn
int bessel_yn(int N, double *z, int nZ, double *y_n, double *dy_n)
Computes the values of the spherical Bessel function of the second kind (yn) and it's derivative (dyn...
Definition saf_utility_bessel.c:818
bessel_jn
int bessel_jn(int N, double *z, int nZ, double *j_n, double *dj_n)
Computes the values of the spherical Bessel function of the first kind (jn) and it's derivative (djn)
Definition saf_utility_bessel.c:630
cart2sph
void cart2sph(float *cart, int nDirs, int anglesInDegreesFLAG, float *sph)
Converts Cartesian coordinates to spherical coordinates.
Definition saf_utility_geometry.c:298
qmf_create
void qmf_create(void **const phQMF, int nCHin, int nCHout, int hopsize, int hybridmode, QMF_FDDATA_FORMAT format)
Creates an instance of the qmf filterbank.
Definition saf_utility_qmf.c:150
bessel_kn
int bessel_kn(int N, double *z, int nZ, double *k_n, double *dk_n)
Computes the values of the modified spherical Bessel function of the second kind (kn) and it's deriva...
Definition saf_utility_bessel.c:907
evalIIRTransferFunction
void evalIIRTransferFunction(double *b_coeff, double *a_coeff, int nCoeffs, float *freqs, int nFreqs, float fs, int mag2dB, float *magnitude, float *phase_rad)
Computes magnitude and phase response of an IIR filter from its coefficients at user-specified freque...
Definition saf_utility_filters.c:674
latticeDecorrelator_destroy
void latticeDecorrelator_destroy(void **phDecor)
Destroys an instance of the lattice all-pass-filter-based multi-channel signal decorrelator.
Definition saf_utility_decor.c:296
L2_norm
float L2_norm(float *v, int lenV)
Returns the L2 (Euclidean) norm of an arbitrary length vector.
Definition saf_utility_geometry.c:423
utility_svvdot
void utility_svvdot(const float *a, const float *b, const int len, float *c)
Single-precision, vector-vector dot product, i.e.
Definition saf_utility_veclib.c:1208
delaunaynd
void delaunaynd(const float *points, const int nPoints, const int nd, int **DT, int *nDT)
Computes the Delaunay triangulation of an arrangement of points in N-dimensional space.
Definition saf_utility_geometry.c:567
bessel_Yn
void bessel_Yn(int N, double *z, int nZ, double *Y_n, double *dY_n)
Computes the values of the (cylindrical) Bessel function of the second kind (Yn) and it's derivative ...
Definition saf_utility_bessel.c:443
saf_stft_flushBuffers
void saf_stft_flushBuffers(void *const hSTFT)
Flushes the internal buffers with zeros.
Definition saf_utility_fft.c:480
saf_fft_destroy
void saf_fft_destroy(void **const phFFT)
Destroys an instance of saf_fft.
Definition saf_utility_fft.c:844
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
saf_stft_create
void saf_stft_create(void **const phSTFT, int winsize, int hopsize, int nCHin, int nCHout, SAF_STFT_FDDATA_FORMAT FDformat)
Creates an instance of saf_stft.
Definition saf_utility_fft.c:298
calcDVFShelfParams
void calcDVFShelfParams(int i, float rhoIn, float *g0, float *gInf, float *fc)
Calculate the high shelf gains and cutoff parameters, given a azimuth index i and distance rho.
Definition saf_utility_dvf.c:82
bessel_in
int bessel_in(int N, double *z, int nZ, double *i_n, double *di_n)
Computes the values of the modified spherical Bessel function of the first kind (in) and it's derivat...
Definition saf_utility_bessel.c:724
faf_IIRFilterbank_apply
void faf_IIRFilterbank_apply(void *hFaF, float *inSig, float **outBands, int nSamples)
Applies the Favrot & Faller filterbank.
Definition saf_utility_filters.c:1176
BUTTER_FILTER_BSF
@ BUTTER_FILTER_BSF
band-stop filter
Definition saf_utility_filters.h:70
BUTTER_FILTER_LPF
@ BUTTER_FILTER_LPF
low-pass filter
Definition saf_utility_filters.h:67
BUTTER_FILTER_HPF
@ BUTTER_FILTER_HPF
high-pass filter
Definition saf_utility_filters.h:68
BUTTER_FILTER_BPF
@ BUTTER_FILTER_BPF
band-pass filter
Definition saf_utility_filters.h:69
EULER_ROTATION_YAW_PITCH_ROLL
@ EULER_ROTATION_YAW_PITCH_ROLL
yaw-pitch-roll, 'zyx'
Definition saf_utility_geometry.h:51
SAF_STFT_BANDS_CH_TIME
@ SAF_STFT_BANDS_CH_TIME
nBands x nChannels x nTimeHops
Definition saf_utility_fft.h:50
QMF_BANDS_CH_TIME
@ QMF_BANDS_CH_TIME
nBands x nChannels x nTimeHops
Definition saf_utility_qmf.h:44
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
calloc1d
void * calloc1d(size_t dim1, size_t data_size)
1-D calloc (same as calloc, but with error checking)
Definition md_malloc.c:69
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.
quaternion_data
Quaternion data structure.
Definition saf_utility_geometry.h:36
quaternion_data::Q
float Q[4]
WXYZ values of the quaternion [-1..1].
Definition saf_utility_geometry.h:43
quaternion_data::x
float x
X value of the quaternion [-1..1].
Definition saf_utility_geometry.h:39
quaternion_data::z
float z
Z value of the quaternion [-1..1].
Definition saf_utility_geometry.h:41
quaternion_data::y
float y
Y value of the quaternion [-1..1].
Definition saf_utility_geometry.h:40
quaternion_data::w
float w
W value of the quaternion [-1..1].
Definition saf_utility_geometry.h:38
test__saf_fft
void test__saf_fft(void)
Testing the forward and backward complex-complex FFT (saf_fft)
Definition test__utilities_module.c:413
test__evalIIRTransferFunction
void test__evalIIRTransferFunction(void)
Testing that the magnitudes and phases returned by evalIIRTransferFunction() are numerically similar ...
Definition test__utilities_module.c:993
test__getVoronoiWeights
void test__getVoronoiWeights(void)
Testing that the weights from the getVoronoiWeights() function sum to 4pi and that the weights are al...
Definition test__utilities_module.c:681
test__sphericalBesselFunctions
void test__sphericalBesselFunctions(void)
Testing bessel_jn(), bessel_in(), bessel_yn(), and bessel_kn()
Definition test__utilities_module.c:54
test__saf_stft_50pc_overlap
void test__saf_stft_50pc_overlap(void)
Testing for perfect reconstruction of the saf_stft (when configured for 50% window overlap)
Definition test__utilities_module.c:207
test__butterCoeffs
void test__butterCoeffs(void)
Testing that the coefficients computed with butterCoeffs() are numerically similar to the "butter" fu...
Definition test__utilities_module.c:870
test__sortf
void test__sortf(void)
Testing the sortf() function (sorting real floating point numbers)
Definition test__utilities_module.c:560
test__cmplxPairUp
void test__cmplxPairUp(void)
Testing the cmplxPairUp() function (grouping up conjugate symmetric values)
Definition test__utilities_module.c:638
test__smb_pitchShifter
void test__smb_pitchShifter(void)
Testing that the smb_pitchShifter can shift the energy of input spectra by one octave down.
Definition test__utilities_module.c:520
test__unique_i
void test__unique_i(void)
Testing that the unique_i() function operates correctly.
Definition test__utilities_module.c:736
test__saf_rfft
void test__saf_rfft(void)
Testing the forward and backward real-(half)complex FFT (saf_rfft)
Definition test__utilities_module.c:374
test__cylindricalBesselFunctions
void test__cylindricalBesselFunctions(void)
Testing bessel_Jn(), bessel_Yn()
Definition test__utilities_module.c:27
test__dvf_dvfShelfCoeffs
void test__dvf_dvfShelfCoeffs(void)
Test the generation of high shelf coeffs based on shelf gain and fc parameters.
Definition test__utilities_module.c:1398
test__delaunaynd
void test__delaunaynd(void)
Testing that the delaunaynd() function can triangulate basic shapes.
Definition test__utilities_module.c:123
test__faf_IIRFilterbank
void test__faf_IIRFilterbank(void)
Testing that the faf_IIRFilterbank can reconstruct the original signal power.
Definition test__utilities_module.c:1195
test__qmf
void test__qmf(void)
Testing the (near)-perfect reconstruction performance of the QMF filterbank.
Definition test__utilities_module.c:454
test__gexpm
void test__gexpm(void)
Testing computing the matrix exponential - comparing the output to that of the "expm" function in Mat...
Definition test__utilities_module.c:1268
test__dvf_interpDVFShelfParams
void test__dvf_interpDVFShelfParams(void)
Test the interpolation of high shelf parameters based on distance and incidence angle parameters.
Definition test__utilities_module.c:1348
test__sortz
void test__sortz(void)
Testing the sortz() function (sorting complex double-floating point numbers)
Definition test__utilities_module.c:592
test__cart2sph
void test__cart2sph(void)
Testing cart2sph() and sph2cart() are reversible.
Definition test__utilities_module.c:99
test__dvf_calcDVFShelfParams
void test__dvf_calcDVFShelfParams(void)
Calculate high shelf parameters, g0, gInf, fc, from the lookup table coefficients (10 degree steps)
Definition test__utilities_module.c:1304
test__saf_matrixConv
void test__saf_matrixConv(void)
Testing the saf_matrixConv.
Definition test__utilities_module.c:330
test__latticeDecorrelator
void test__latticeDecorrelator(void)
Testing the performance of the latticeDecorrelator, verifying that the inter- channel cross-correlati...
Definition test__utilities_module.c:768
test__saf_stft_LTI
void test__saf_stft_LTI(void)
Testing for perfect reconstruction of the saf_stft (when configured for linear time-invariant (LTI) f...
Definition test__utilities_module.c:271
test__quaternion
void test__quaternion(void)
Testing that quaternion2rotationMatrix() and rotationMatrix2quaternion() are reversible.
Definition test__utilities_module.c:170