Looking for MKL FFT complementary functions

[TripleS]

Hey 

I migrated my DFT code from IPP to FFT code using MKL FFT functions, I'm looking forward for performance improvement.

1st, I want to do phase correlation in spectrum domain. 

In my algorithm I have 2 input images and the output is a correlation surface 

I do the following: 

1. FFT both input images 

2. for each element in the spatial domain, I'd like to multiply with complex conjugate, than normalize  each element: 

2.1       (a1 + ib1) * (a2 - ib2) 

2.2       (a + ib) / magnidute(a,b)  --> Z 

3. inverse FFT Z 

I'm looking for complementary functions for FFT library that calculate a multiplication with complex conjugate, and function that calculate each element magnitude.  

I thought I'd be great to use: 

vcMulByConj 

and 

VCABS 

functions that reside on "mkl_vml_functions.h" 

However, here comes a problem of using MKL FFT memory model (DFTI_CONJUGATE_EVEN_STORAGE  in packed memory format ), I don't know how to access each element in the spatial domain, how to build the matrix back in packed format, and most important is how to make this thing work extremely fast.  

I'm looking for functions that already know how to work with packed format and execute such functionality on packed format matrix, 

Are there functions that gives such functionality? 

Do I need to write this functionally by myself ? How to access each matrix element ? 

2nd, I'm also want to do FFT shift in the spatial domain by doing (-1)^(i+j) for each element. Is there a function that can do such thing (perhaps matrix pixel wise multiplication) ? Working with packed format. 

Thanks 

S 

[Ying_H_Intel]

Hi 

I guess, you are doing some convolution operation on Image, right?  The performance may depend on many factors, like your image size etc.

1) About packed format

If you do it with MKL FFT API.  

for example, if with DftiSetValue(h1,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX); you may don't need to take care of the packed model. 

If for the package format, you can refer to MKL manual in MKL documention (https://software.intel.com/en-us/articles/intel-math-kernel-library-documentation)  => Section

DFTI_CCS_FORMAT for Two-dimensional Transforms

The following figure illustrates the storage of a two-dimensional (2D) M-by-N conjugate-even sequence in a
real array for the CCS packed format. This format requires an array of size (M+2)-by-(N+2). Row-major
layout and zero-based indexing are used. Different colors mark logically separate parts of the result. "n/u"
means "not used".
Storage of a 2D M-by-N Conjugate-even Sequence in a Real Array for the CCS Format

2) yes, you can use the function to do (a + ib) / magnidute(a,b)

vcMulByConj 

and 

VCABS  or v?LinearFrac

3)  FFT shift in the spatial domain by doing (-1)^(i+j) for each element. Is there a function that can do such thing (perhaps matrix pixel wise multiplication) ?

yes, maybe matrix pixel wise multiplication by VML function,  or consider FFT DFTI_FORWARD_SCALE, which allow you multiply some value to spatial. 

Here is one use FFT for image reconstruction sample, for your reference. 

https://software.intel.com/en-us/node/507042

Best Regards,

Ying

[TripleS]

Hey Ying 

Thank you for your reply 

I'd like to use: DftiSetValue(h1,DFTI_CONJUGATE_EVEN_STORAGE,DFTI_COMPLEX_COMPLEX)  which don't take care of the packed memory format. 

I expected that the memory would save in conjugate even storage, which means for the number of elements in read and spatial domain the following: 

for Input size: sizeof(float) * M * N  

The output size would be : sizeof( complex float ) * M * ( N / 2 + 1 ) 

Which is exactly what am I trying to do, Please have a look at my following example, I'm trying to figure out what am I missing here. 

DftiComputeForward throws an exception. 

( Could be an output stride ? I thought the ouput stride should set by default to default parameters ) 

Thanks 

#include <math.h>
#include <mkl.h>
#include <mkl_dfti.h>
#include <mkl_cblas.h>

using namespace std;


static int Test()
{
	int m_nWorkingRows = 1024;
	int m_nWorkingCols = 1024;

	MKL_LONG mklStatus;
	DFTI_DESCRIPTOR_HANDLE handleDFTI_foward;

	MKL_LONG descSize[2] = { m_nWorkingRows, m_nWorkingCols };

	// Note: the output is DFTI_EVEN_STORAGE. store only the conjugate even domain by default, real data scope is d2 * d1, imaginary data scope is: d2 * (d1 / 2 + 1)
	mklStatus = DftiCreateDescriptor(&handleDFTI_foward, DFTI_SINGLE, DFTI_REAL, 2, descSize);
	if (mklStatus && !DftiErrorClass(mklStatus, DFTI_NO_ERROR)) { /*printf("Error: %s\n", DftiErrorMessage(mklStatus));*/  return false; }

	mklStatus = DftiSetValue(handleDFTI_foward, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
	if (mklStatus && !DftiErrorClass(mklStatus, DFTI_NO_ERROR)) { /*printf("Error: %s\n", DftiErrorMessage(mklStatus));*/  return false; }

	mklStatus = DftiSetValue(handleDFTI_foward, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
	if (mklStatus && !DftiErrorClass(mklStatus, DFTI_NO_ERROR)) { /*printf("Error: %s\n", DftiErrorMessage(mklStatus));*/  return false; }

	mklStatus = DftiCommitDescriptor(handleDFTI_foward);
	if (mklStatus && !DftiErrorClass(mklStatus, DFTI_NO_ERROR)) { /*printf("Error: %s\n", DftiErrorMessage(mklStatus));*/  return false; }


	
	//// 
	void* srcData = MKL_malloc(sizeof(float) * m_nWorkingRows * m_nWorkingCols, 64);
	void* dstData = MKL_malloc(sizeof(float) * 2 * m_nWorkingRows * (m_nWorkingCols / 2 + 1), 64);


	// THE FOLLOWING LINE THROWS AN EXCEPTION: Unhandled exception at 0x000007FEC9B7233C (mkl_mc3.dll) in xxxxx.exe: 0xC0000005: Access violation writing location 0x000000000B7BE080.
	mklStatus = DftiComputeForward(handleDFTI_foward, srcData, dstData);
	if (mklStatus && !DftiErrorClass(mklStatus, DFTI_NO_ERROR)) { /*printf("Error: %s\n", DftiErrorMessage(mklStatus));*/  return false; }

	return true;
}