Function nc_crosspower

numexpr3/complex_functions.hpp:1672–1703

RAM: cross-power-spectrum (for Phase Correlation). https://en.wikipedia.org/wiki/Phase_correlation

// RAM: cross-power-spectrum (for Phase Correlation).
// https://en.wikipedia.org/wiki/Phase_correlation
static void
nc_crosspower(npy_intp n, npy_complex64 *a, npy_intp sb1, npy_complex64 *b, npy_intp sb2, npy_complex64 *r) {
    npy_float32 mag;
    npy_complex64 prod; // Not-inlining this into r should improve SIMD performance.

    if( sb1 == sizeof(npy_complex64) && sb2 == sizeof(npy_complex64) ) { // Aligned
        for( npy_intp I = 0; I < n; I++ ) {
            // Take product of A and complex conjugate of B
            prod.real =  a[I].real*b[I].real + a[I].imag*b[I].imag;
            prod.imag = -a[I].real*b[I].imag + a[I].imag*b[I].real;
            // Compute magnitude 
            mag = 1.0f/sqrtf(prod.real*prod.real + prod.imag*prod.imag);
            // Normalize
            r[I].real = prod.real*mag;
            r[I].imag = prod.imag*mag;
        }
    }
    else {
        sb1 /= sizeof(npy_complex64);
        sb2 /= sizeof(npy_complex64);
        for( npy_intp I = 0; I < n; I++ ) {
            // Take product of A and complex conjugate of B
            prod.real =  a[I*sb1].real*b[I*sb2].real + a[I*sb1].imag*b[I*sb2].imag;
            prod.imag = -a[I*sb1].real*b[I*sb2].imag + a[I*sb1].imag*b[I*sb2].real;
            // Compute magnitude 
            mag = 1.0f/sqrtf(prod.real*prod.real + prod.imag*prod.imag);
            // Normalize
            r[I].real = prod.real*mag;
            r[I].imag = prod.imag*mag;
        }
    }
}

static void
nc_crosspower(npy_intp n, npy_complex128 *a, npy_intp sb1, npy_complex128 *b, npy_intp sb2, npy_complex128 *r) {