SAF
Loading...
Searching...
No Matches
saf_cdf4sap.c
Go to the documentation of this file.
1/*
2 * Copyright 2018-2019 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
40#include "saf_cdf4sap.h"
41#include "saf_externals.h"
42#include "../saf_utilities/saf_utilities.h"
43
48typedef struct _cdf4sap_data {
49 /* Dimensions of Cx and Cy */
50 int nXcols, nYcols;
51
52 /* intermediate vectors & matrices */
53 void* hSVD;
54 float* lambda, *U_Cy, *S_Cy, *Ky, *U_Cx, *S_Cx, *s_Cx, *Kx, *Kx_reg_inverse, *U, *V, *P;
55 float* G_hat, *Cx_QH;
56 float* GhatH_Ky, *QH_GhatH_Ky, *KxH_QH_GhatH_Ky, *lambda_UH;
57 float* P_Kxreginverse;
58 float* Cx_MH, *Cy_tilde;
59 float* G_M;
60
61}cdf4sap_data;
62
67typedef struct _cdf4sap_cmplx_data {
68 /* Dimensions of Cx and Cy */
69 int nXcols, nYcols;
70
71 /* intermediate vectors & matrices */
72 void* hSVD;
73 float_complex* Cr_cmplx;
74 float_complex* lambda, *U_Cy, *S_Cy, *S_Cx, *Ky, *U_Cx, *Kx, *Kx_reg_inverse, *U, *V, *P;
75 float* s_Cx, *G_hat_diag;
76 float_complex* G_hat, *Cx_QH;
77 float_complex* GhatH_Ky, *QH_GhatH_Ky, *KxH_QH_GhatH_Ky, *lambda_UH;
78 float_complex *P_Kxreginverse;
79 float_complex *Cx_MH, *Cy_tilde;
80 float_complex* G_M;
81
82}cdf4sap_cmplx_data;
83
84void cdf4sap_create
85(
86 void ** const phCdf,
87 int nXcols,
88 int nYcols
89)
90{
91 *phCdf = malloc1d(sizeof(cdf4sap_data));
92 cdf4sap_data *h = (cdf4sap_data*)(*phCdf);
93
94 h->nXcols = nXcols;
95 h->nYcols = nYcols;
96 h->lambda = malloc1d(nYcols * nXcols * sizeof(float));
97
98 /* For the SVD */
99 utility_ssvd_create(&h->hSVD, SAF_MAX(nXcols, nYcols), SAF_MAX(nXcols, nYcols));
100
101 /* For the decomposition of Cy */
102 h->U_Cy = malloc1d(nYcols*nYcols*sizeof(float));
103 h->S_Cy = malloc1d(nYcols*nYcols*sizeof(float));
104 h->Ky = malloc1d(nYcols*nYcols*sizeof(float));
105
106 /* For the decomposition of Cx */
107 h->U_Cx = malloc1d(nXcols*nXcols*sizeof(float));
108 h->S_Cx = malloc1d(nXcols*nXcols*sizeof(float));
109 h->s_Cx = malloc1d(nXcols*sizeof(float));
110 h->Kx = malloc1d(nXcols*nXcols*sizeof(float));
111
112 /* For the formulation of regularised Kx^-1 */
113 h->Kx_reg_inverse = malloc1d(nXcols*nXcols*sizeof(float));
114
115 /* For the formulation of normalisation matrix G_hat */
116 h->G_hat = malloc1d(nYcols*nYcols*sizeof(float));
117 h->Cx_QH = malloc1d(nXcols*nYcols*sizeof(float));
118
119 /* For the formulation of optimal P */
120 h->GhatH_Ky = malloc1d(nYcols*nYcols*sizeof(float));
121 h->QH_GhatH_Ky = malloc1d(nXcols*nYcols*sizeof(float));
122 h->KxH_QH_GhatH_Ky = malloc1d(nXcols*nYcols*sizeof(float));
123 h->U = malloc1d(nXcols*nXcols*sizeof(float));
124 h->V = malloc1d(nYcols*nYcols*sizeof(float));
125 h->lambda_UH = malloc1d(nYcols*nXcols*sizeof(float));
126 h->P = malloc1d(nYcols*nXcols*sizeof(float));
127
128 /* For the formulation of M */
129 h->P_Kxreginverse = malloc1d(nYcols*nXcols*sizeof(float));
130
131 /* For the formulation of the residual covariance matrix */
132 h->Cx_MH = malloc1d(nXcols*nYcols*sizeof(float));
133 h->Cy_tilde = malloc1d(nYcols*nYcols*sizeof(float));
134
135 /* For using energy compensation instead of residuals */
136 h->G_M = malloc1d(nYcols*nXcols*sizeof(float));
137}
138
139void cdf4sap_cmplx_create
140(
141 void ** const phCdf,
142 int nXcols,
143 int nYcols
144)
145{
146 *phCdf = malloc1d(sizeof(cdf4sap_cmplx_data));
147 cdf4sap_cmplx_data *h = (cdf4sap_cmplx_data*)(*phCdf);
148
149 h->nXcols = nXcols;
150 h->nYcols = nYcols;
151 h->lambda = malloc1d(nYcols * nXcols * sizeof(float_complex));
152 h->Cr_cmplx = malloc1d(nYcols * nYcols * sizeof(float_complex));
153
154 /* For the SVD */
155 utility_csvd_create(&h->hSVD, SAF_MAX(nXcols, nYcols), SAF_MAX(nXcols, nYcols));
156
157 /* For the decomposition of Cy */
158 h->U_Cy = malloc1d(nYcols*nYcols*sizeof(float_complex));
159 h->S_Cy = malloc1d(nYcols*nYcols*sizeof(float_complex));
160 h->Ky = malloc1d(nYcols*nYcols*sizeof(float_complex));
161
162 /* For the decomposition of Cx */
163 h->U_Cx = malloc1d(nXcols*nXcols*sizeof(float_complex));
164 h->S_Cx = malloc1d(nXcols*nXcols*sizeof(float_complex));
165 h->s_Cx = malloc1d(nXcols*sizeof(float));
166 h->Kx = malloc1d(nXcols*nXcols*sizeof(float_complex));
167
168 /* For the formulation of regularised Kx^-1 */
169 h->Kx_reg_inverse = malloc1d(nXcols*nXcols*sizeof(float_complex));
170
171 /* For the formulation of normalisation matrix G_hat */
172 h->G_hat_diag = malloc1d(nYcols*sizeof(float));
173 h->G_hat = malloc1d(nYcols*nYcols*sizeof(float_complex));
174 h->Cx_QH = malloc1d(nXcols*nYcols*sizeof(float_complex));
175
176 /* For the formulation of optimal P */
177 h->GhatH_Ky = malloc1d(nYcols*nYcols*sizeof(float_complex));
178 h->QH_GhatH_Ky = malloc1d(nXcols*nYcols*sizeof(float_complex));
179 h->KxH_QH_GhatH_Ky = malloc1d(nXcols*nYcols*sizeof(float_complex));
180 h->U = malloc1d(nXcols*nXcols*sizeof(float_complex));
181 h->V = malloc1d(nYcols*nYcols*sizeof(float_complex));
182 h->lambda_UH = malloc1d(nYcols*nXcols*sizeof(float_complex));
183 h->P = malloc1d(nYcols*nXcols*sizeof(float_complex));
184
185 /* For the formulation of M */
186 h->P_Kxreginverse = malloc1d(nYcols*nXcols*sizeof(float_complex));
187
188 /* For the formulation of the residual covariance matrix */
189 h->Cx_MH = malloc1d(nXcols*nYcols*sizeof(float_complex));
190 h->Cy_tilde = malloc1d(nYcols*nYcols*sizeof(float_complex));
191
192 /* For using energy compensation instead of residuals */
193 h->G_M = malloc1d(nYcols*nXcols*sizeof(float_complex));
194}
195
196void cdf4sap_destroy
197(
198 void ** const phCdf
199)
200{
201 cdf4sap_data *h = (cdf4sap_data*)(*phCdf);
202
203 if(h!=NULL){
204 utility_ssvd_destroy(&h->hSVD);
205 free(h->lambda);
206 free(h->U_Cy);
207 free(h->S_Cy);
208 free(h->Ky);
209 free(h->U_Cx);
210 free(h->S_Cx);
211 free(h->s_Cx);
212 free(h->Kx);
213 free(h->Kx_reg_inverse);
214 free(h->G_hat);
215 free(h->Cx_QH);
216 free(h->GhatH_Ky);
217 free(h->QH_GhatH_Ky);
218 free(h->KxH_QH_GhatH_Ky);
219 free(h->U);
220 free(h->V);
221 free(h->lambda_UH);
222 free(h->P);
223 free(h->P_Kxreginverse);
224 free(h->Cx_MH);
225 free(h->Cy_tilde);
226 free(h->G_M);
227 free(h);
228 h = NULL;
229 *phCdf = NULL;
230 }
231}
232
233void cdf4sap_cmplx_destroy
234(
235 void ** const phCdf
236)
237{
238 cdf4sap_cmplx_data *h = (cdf4sap_cmplx_data*)(*phCdf);
239
240 if(h!=NULL){
241 utility_csvd_destroy(&h->hSVD);
242 free(h->lambda);
243 free(h->Cr_cmplx);
244 free(h->U_Cy);
245 free(h->S_Cy);
246 free(h->Ky);
247 free(h->U_Cx);
248 free(h->S_Cx);
249 free(h->s_Cx);
250 free(h->Kx);
251 free(h->Kx_reg_inverse);
252 free(h->G_hat_diag);
253 free(h->G_hat);
254 free(h->Cx_QH);
255 free(h->GhatH_Ky);
256 free(h->QH_GhatH_Ky);
257 free(h->KxH_QH_GhatH_Ky);
258 free(h->U);
259 free(h->V);
260 free(h->lambda_UH);
261 free(h->P);
262 free(h->P_Kxreginverse);
263 free(h->Cx_MH);
264 free(h->Cy_tilde);
265 free(h->G_M);
266 free(h);
267 h = NULL;
268 *phCdf = NULL;
269 }
270}
271
272void formulate_M_and_Cr
273(
274 void * const hCdf,
275 float* Cx,
276 float* Cy,
277 float* Q,
278 int useEnergyFLAG,
279 float reg,
280 float* M,
281 float* Cr
282)
283{
284 cdf4sap_data *h = (cdf4sap_data*)(hCdf);
285 int i, j, nXcols, nYcols;
286
287 nXcols = h->nXcols;
288 nYcols = h->nYcols;
289
290 memset(h->lambda, 0, nYcols * nXcols * sizeof(float));
291 for(i = 0; i<SAF_MIN(nXcols,nYcols); i++)
292 h->lambda[i*nXcols + i] = 1.0f;
293
294 /* Decomposition of Cy */
295 utility_ssvd(h->hSVD, Cy, nYcols, nYcols, h->U_Cy, h->S_Cy, NULL, NULL);
296 for(i=0; i< nYcols; i++)
297 h->S_Cy[i*nYcols+i] = sqrtf(SAF_MAX(h->S_Cy[i*nYcols+i], 2.23e-20f));
298 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nYcols, nYcols, 1.0f,
299 h->U_Cy, nYcols,
300 h->S_Cy, nYcols, 0.0f,
301 h->Ky, nYcols);
302
303 /* Decomposition of Cx */
304 utility_ssvd(h->hSVD, Cx, nXcols, nXcols, h->U_Cx, h->S_Cx, NULL, h->s_Cx);
305 for(i=0; i< nXcols; i++){
306 h->S_Cx[i*nXcols+i] = sqrtf(SAF_MAX(h->S_Cx[i*nXcols+i], 2.23e-20f));
307 h->s_Cx[i] = sqrtf(SAF_MAX(h->s_Cx[i], 2.23e-20f));
308 }
309 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nXcols, nXcols, nXcols, 1.0f,
310 h->U_Cx, nXcols,
311 h->S_Cx, nXcols, 0.0f,
312 h->Kx, nXcols);
313
314 /* Regularisation of S_Cx */
315 int ind;
316 float limit, maxVal;
317 utility_simaxv(h->s_Cx, nXcols, &ind);
318 limit = h->s_Cx[ind] * reg + 2.23e-13f;
319 for(i=0; i < nXcols; i++)
320 h->S_Cx[i*nXcols+i] = 1.0f / SAF_MAX(h->S_Cx[i*nXcols+i], limit);
321
322 /* Formulate regularised Kx^-1 */
323 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, nXcols, nXcols, nXcols, 1.0f,
324 h->S_Cx, nXcols,
325 h->U_Cx, nXcols, 0.0f,
326 h->Kx_reg_inverse, nXcols);
327
328 /* Formulate normalisation matrix G_hat */
329 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, nXcols, nYcols, nXcols, 1.0f,
330 Cx, nXcols,
331 Q, nXcols, 0.0f,
332 h->Cx_QH, nYcols);
333 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nYcols, nXcols, 1.0f,
334 Q, nXcols,
335 h->Cx_QH, nYcols, 0.0f,
336 h->G_hat, nYcols);
337 maxVal = -2.23e13f;
338 for(i=0; i< nYcols; i++)
339 maxVal = h->G_hat[i*nYcols+i] > maxVal ? h->G_hat[i*nYcols+i] : maxVal;
340 limit = maxVal * 0.001f + 2.23e-13f;
341 for(i=0; i < nYcols; i++)
342 for(j=0; j < nYcols; j++)
343 h->G_hat[i*nYcols+j] = i==j ? sqrtf(SAF_MAX(Cy[i*nYcols+j], 2.23e-13f) / SAF_MAX(h->G_hat[i*nYcols+j], limit)) : 0.0f;
344
345 /* Formulate optimal P */
346 cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, nYcols, nYcols, nYcols, 1.0f,
347 h->G_hat, nYcols,
348 h->Ky, nYcols, 0.0f,
349 h->GhatH_Ky, nYcols);
350 cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, nXcols, nYcols, nYcols, 1.0f,
351 Q, nXcols,
352 h->GhatH_Ky, nYcols, 0.0f,
353 h->QH_GhatH_Ky, nYcols);
354 cblas_sgemm(CblasRowMajor, CblasTrans, CblasNoTrans, nXcols, nYcols, nXcols, 1.0f,
355 h->Kx, nXcols,
356 h->QH_GhatH_Ky, nYcols, 0.0f,
357 h->KxH_QH_GhatH_Ky, nYcols);
358 utility_ssvd(h->hSVD, h->KxH_QH_GhatH_Ky, nXcols, nYcols, h->U, NULL, h->V, NULL);
359 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, nYcols, nXcols, nXcols, 1.0f,
360 h->lambda, nXcols,
361 h->U, nXcols, 0.0f,
362 h->lambda_UH, nXcols);
363 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nXcols, nYcols, 1.0f,
364 h->V, nYcols,
365 h->lambda_UH, nXcols, 0.0f,
366 h->P, nXcols);
367
368 /* Formulate M */
369 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nXcols, nXcols, 1.0f,
370 h->P, nXcols,
371 h->Kx_reg_inverse, nXcols, 0.0f,
372 h->P_Kxreginverse, nXcols);
373 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nXcols, nYcols, 1.0f,
374 h->Ky, nYcols,
375 h->P_Kxreginverse, nXcols, 0.0f,
376 M, nXcols);
377
378 /* Formulate residual covariance matrix */
379 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, nXcols, nYcols, nXcols, 1.0f,
380 Cx, nXcols,
381 M, nXcols, 0.0f,
382 h->Cx_MH, nYcols);
383 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nYcols, nXcols, 1.0f,
384 M, nXcols,
385 h->Cx_MH, nYcols, 0.0f,
386 h->Cy_tilde, nYcols);
387 if(Cr != NULL)
388 for(i=0; i < nYcols*nYcols; i++)
389 Cr[i] = Cy[i] - h->Cy_tilde[i];
390
391 /* Use energy compensation instead of residuals */
392 if(useEnergyFLAG){
393 for(i=0; i < nYcols; i++)
394 for(j=0; j < nYcols; j++)
395 h->G_hat[i*nYcols+j] = i==j ? sqrtf( SAF_MAX(Cy[i*nYcols+j], 2.23e-20f) / (h->Cy_tilde[i*nYcols+j]+2.23e-7f)) : 0.0f;
396 cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nXcols, nYcols, 1.0f,
397 h->G_hat, nYcols,
398 M, nXcols, 0.0f,
399 h->G_M, nXcols);
400 memcpy(M, h->G_M, nYcols*nXcols*sizeof(float));
401 if(Cr != NULL)
402 memset(Cr, 0, nYcols*nYcols*sizeof(float));
403 }
404}
405
406void formulate_M_and_Cr_cmplx
407(
408 void * const hCdf,
409 float_complex* Cx,
410 float_complex* Cy,
411 float_complex* Q,
412 int useEnergyFLAG,
413 float reg,
414 float_complex* M,
415 float_complex* Cr
416)
417{
418 cdf4sap_cmplx_data *h = (cdf4sap_cmplx_data*)(hCdf);
419 int i, j, nXcols, nYcols;
420 const float_complex calpha = cmplxf(1.0f, 0.0f); const float_complex cbeta = cmplxf(0.0f, 0.0f);
421
422 nXcols = h->nXcols;
423 nYcols = h->nYcols;
424
425 memset(h->lambda, 0, nYcols * nXcols * sizeof(float_complex));
426 for(i = 0; i<SAF_MIN(nXcols,nYcols); i++)
427 h->lambda[i*nXcols + i] = cmplxf(1.0f, 0.0f);
428
429 /* Decomposition of Cy */
430 utility_csvd(h->hSVD, Cy, nYcols, nYcols, h->U_Cy, h->S_Cy, NULL, NULL);
431 for(i=0; i< nYcols; i++)
432 h->S_Cy[i*nYcols+i] = cmplxf(sqrtf(SAF_MAX(crealf(h->S_Cy[i*nYcols+i]), 2.23e-20f)), 0.0f);
433 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nYcols, nYcols, &calpha,
434 h->U_Cy, nYcols,
435 h->S_Cy, nYcols, &cbeta,
436 h->Ky, nYcols);
437
438 /* Decomposition of Cx */
439 utility_csvd(h->hSVD, Cx, nXcols, nXcols, h->U_Cx, h->S_Cx, NULL, h->s_Cx);
440 for(i=0; i< nXcols; i++){
441 h->s_Cx[i] = sqrtf(SAF_MAX(h->s_Cx[i], 2.23e-13f));
442 h->S_Cx[i*nXcols+i] = cmplxf(h->s_Cx[i], 0.0f);
443 }
444 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nXcols, nXcols, nXcols, &calpha,
445 h->U_Cx, nXcols,
446 h->S_Cx, nXcols, &cbeta,
447 h->Kx, nXcols);
448
449 /* Regularisation of S_Cx */
450 int ind;
451 float limit, maxVal;
452 //utility_simaxv(h->s_Cx, nXcols, &ind);
453 ind = 0; /* utility_csvd returns the singular values in decending order */
454 limit = h->s_Cx[ind] * reg + 2.23e-13f;
455 for(i=0; i < nXcols; i++)
456 h->S_Cx[i*nXcols+i] = cmplxf(1.0f / SAF_MAX(h->s_Cx[i], limit), 0.0f);
457
458 /* Formulate regularised Kx^-1 */
459 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, nXcols, nXcols, nXcols, &calpha,
460 h->S_Cx, nXcols,
461 h->U_Cx, nXcols, &cbeta,
462 h->Kx_reg_inverse, nXcols);
463
464 /* Formulate normalisation matrix G_hat */
465 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, nXcols, nYcols, nXcols, &calpha,
466 Cx, nXcols,
467 Q, nXcols, &cbeta,
468 h->Cx_QH, nYcols);
469 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nYcols, nXcols, &calpha,
470 Q, nXcols,
471 h->Cx_QH, nYcols, &cbeta,
472 h->G_hat, nYcols);
473 /* h->G_hat: imaginary parts along the diagonal are ~0, so it's OK to take
474 * the real instead below: */
475 maxVal = -2.23e13f;
476 for(i=0; i< nYcols; i++)
477 maxVal = cabsf(h->G_hat[i*nYcols+i]) > maxVal ? cabsf(h->G_hat[i*nYcols+i]) : maxVal; // crealf->cabsf
478 limit = maxVal * 0.001f + 2.23e-13f;
479#if 0 /* DOES NOT PASS UNIT TEST: */
480 cblas_scopy(nYcols, (float*)h->G_hat, 2*(nYcols+1), h->G_hat_diag, 1); /* take diagonal (real) */
481 memset(h->G_hat, 0, nYcols*nYcols*sizeof(float_complex));
482 for(i=0; i<nYcols; i++)
483 h->G_hat_diag[i] = sqrtf(((float*)Cy)[2*(i*nYcols+i)]/SAF_MAX(h->G_hat_diag[i], limit));
484 cblas_scopy(nYcols, (float*)h->G_hat_diag, 1, (float*)h->G_hat, /*re+im*/2*(nYcols+1)); /* load the diagonal */
485#else
486 for(i=0; i < nYcols; i++)
487 for(j=0; j < nYcols; j++)
488 h->G_hat[i*nYcols+j] = i==j ? cmplxf(crealf(csqrtf( ccdivf(Cy[i*nYcols+j], cmplxf(SAF_MAX(cabsf(h->G_hat[i*nYcols+j]), limit), 0.0f)))), 0.0f) : cmplxf(0.0f, 0.0f); // changed crealf->cabsf
489#endif
490
491 /* Formulate optimal P */
492 cblas_cgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans, nYcols, nYcols, nYcols, &calpha,
493 h->G_hat, nYcols,
494 h->Ky, nYcols, &cbeta,
495 h->GhatH_Ky, nYcols);
496 cblas_cgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans, nXcols, nYcols, nYcols, &calpha,
497 Q, nXcols,
498 h->GhatH_Ky, nYcols, &cbeta,
499 h->QH_GhatH_Ky, nYcols);
500 cblas_cgemm(CblasRowMajor, CblasConjTrans, CblasNoTrans, nXcols, nYcols, nXcols, &calpha,
501 h->Kx, nXcols,
502 h->QH_GhatH_Ky, nYcols, &cbeta,
503 h->KxH_QH_GhatH_Ky, nYcols);
504 utility_csvd(h->hSVD, h->KxH_QH_GhatH_Ky, nXcols, nYcols, h->U, NULL, h->V, NULL);
505 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, nYcols, nXcols, nXcols, &calpha,
506 h->lambda, nXcols,
507 h->U, nXcols, &cbeta,
508 h->lambda_UH, nXcols);
509 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nXcols, nYcols, &calpha,
510 h->V, nYcols,
511 h->lambda_UH, nXcols, &cbeta,
512 h->P, nXcols);
513
514 /* Formulate M */
515 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nXcols, nXcols, &calpha,
516 h->P, nXcols,
517 h->Kx_reg_inverse, nXcols, &cbeta,
518 h->P_Kxreginverse, nXcols);
519 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nXcols, nYcols, &calpha,
520 h->Ky, nYcols,
521 h->P_Kxreginverse, nXcols, &cbeta,
522 M, nXcols);
523
524 /* Formulate residual covariance matrix */
525 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasConjTrans, nXcols, nYcols, nXcols, &calpha,
526 Cx, nXcols,
527 M, nXcols, &cbeta,
528 h->Cx_MH, nYcols);
529 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nYcols, nXcols, &calpha,
530 M, nXcols,
531 h->Cx_MH, nYcols, &cbeta,
532 h->Cy_tilde, nYcols);
533 if(Cr != NULL) {
534#if 1
535 cblas_sscal(nYcols*nYcols, 0.0f, ((float*)Cr)+1, 2); /* set imag part to zero */
536 cblas_scopy(nYcols*nYcols, (float*)Cy, 2, (float*)Cr, 2); /* copy real part */
537 cblas_saxpy(nYcols*nYcols, -1.0f, (float*)h->Cy_tilde, 2, (float*)Cr, 2); /* subtract tilde to get residual (real only) */
538#else
539 for(i=0; i < nYcols*nYcols; i++){
540 h->Cr_cmplx[i] = ccsubf(Cy[i], h->Cy_tilde[i]);
541 Cr[i] = cmplxf(crealf(h->Cr_cmplx[i]), 0.0f);
542 //Cr[i] = h->Cr_cmplx[i];
543 }
544#endif
545 }
546
547 /* Use energy compensation instead of residuals */
548 if(useEnergyFLAG){
549 for(i=0; i < nYcols; i++)
550 for(j=0; j < nYcols; j++)
551 h->G_hat[i*nYcols+j] = i==j ? csqrtf(ccdivf(Cy[i*nYcols+j], craddf(h->Cy_tilde[i*nYcols+j], 2.23e-13f))): cmplxf(0.0f, 0.0f);
552
553 /* for(i=0; i < nYcols; i++)
554 for(j=0; j < nYcols; j++)
555 h->G_hat[i*nYcols+j] = i==j ? cmplxf(crealf(csqrtf(ccdivf(Cy[i*nYcols+j], cmplxf(MAX(crealf(h->Cy_tilde[i*nYcols+j])+2.23e-7f, limit), 0.0f)))), 0.0f) : cmplxf(0.0f, 0.0f);
556 */
557 cblas_cgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, nYcols, nXcols, nYcols, &calpha,
558 h->G_hat, nYcols,
559 M, nXcols, &cbeta,
560 h->G_M, nXcols);
561 memcpy(M, h->G_M, nYcols*nXcols*sizeof(float_complex));
562 if(Cr != NULL)
563 memset(Cr, 0, nYcols*nYcols*sizeof(float_complex));
564 }
565}
formulate_M_and_Cr_cmplx
void formulate_M_and_Cr_cmplx(void *const hCdf, float_complex *Cx, float_complex *Cy, float_complex *Q, int useEnergyFLAG, float reg, float_complex *M, float_complex *Cr)
Computes the optimal mixing matrices.
Definition saf_cdf4sap.c:407
cdf4sap_create
void cdf4sap_create(void **const phCdf, int nXcols, int nYcols)
Creates an instance of the Covariance Domain Framework.
Definition saf_cdf4sap.c:85
cdf4sap_cmplx_destroy
void cdf4sap_cmplx_destroy(void **const phCdf)
Destroys an instance of the Covariance Domain Framework.
Definition saf_cdf4sap.c:234
cdf4sap_destroy
void cdf4sap_destroy(void **const phCdf)
Destroys an instance of the Covariance Domain Framework.
Definition saf_cdf4sap.c:197
cdf4sap_cmplx_create
void cdf4sap_cmplx_create(void **const phCdf, int nXcols, int nYcols)
Creates an instance of the Covariance Domain Framework.
Definition saf_cdf4sap.c:140
formulate_M_and_Cr
void formulate_M_and_Cr(void *const hCdf, float *Cx, float *Cy, float *Q, int useEnergyFLAG, float reg, float *M, float *Cr)
Computes the optimal mixing matrices.
Definition saf_cdf4sap.c:273
utility_csvd_destroy
void utility_csvd_destroy(void **const phWork)
De-allocate the working struct used by utility_csvd()
Definition saf_utility_veclib.c:1773
SAF_MAX
#define SAF_MAX(a, b)
Returns the maximum of the two values.
Definition saf_utilities.h:58
utility_ssvd_create
void utility_ssvd_create(void **const phWork, int maxDim1, int maxDim2)
(Optional) Pre-allocate the working struct used by utility_ssvd()
Definition saf_utility_veclib.c:1612
utility_ssvd
void utility_ssvd(void *const hWork, const float *A, const int dim1, const int dim2, float *U, float *S, float *V, float *sing)
Singular value decomposition: single precision, i.e.
Definition saf_utility_veclib.c:1645
utility_csvd
void utility_csvd(void *const hWork, const float_complex *A, const int dim1, const int dim2, float_complex *U, float_complex *S, float_complex *V, float *sing)
Singular value decomposition: single precision complex, i.e.
Definition saf_utility_veclib.c:1792
SAF_MIN
#define SAF_MIN(a, b)
Returns the minimum of the two values.
Definition saf_utilities.h:55
utility_csvd_create
void utility_csvd_create(void **const phWork, int maxDim1, int maxDim2)
(Optional) Pre-allocate the working struct used by utility_csvd()
Definition saf_utility_veclib.c:1757
utility_simaxv
void utility_simaxv(const float *a, const int len, int *index)
Single-precision, index of maximum absolute value in a vector, i.e.
Definition saf_utility_veclib.c:391
utility_ssvd_destroy
void utility_ssvd_destroy(void **const phWork)
De-allocate the working struct used by utility_ssvd()
Definition saf_utility_veclib.c:1627
malloc1d
void * malloc1d(size_t dim1_data_size)
1-D malloc (same as malloc, but with error checking)
Definition md_malloc.c:59
saf_cdf4sap.h
Main header for the Covariance Domain Framework module (SAF_CDF4SAP_MODULE)
saf_externals.h
Include header for SAF externals.
saf_utilities.h
Main header for the utilities module (SAF_UTILITIES_MODULE)
cdf4sap_cmplx_data
Main data structure for the Covariance Domain Framework for Spatial Audio Processing (CDF4SAP),...
Definition saf_cdf4sap.c:67
cdf4sap_data
Main data structure for the Covariance Domain Framework for Spatial Audio Processing (CDF4SAP),...
Definition saf_cdf4sap.c:48