Search code examples
arraysmatlabtypescastingcomplex-numbers

Casting complex to real without data copy in MATLAB R2018a and newer


Since MATLAB R2018a, complex-valued matrices are stored internally as a single data block, with the real and imaginary component of each matrix element stored next to each other -- they call this "interleaved complex". (Previously such matrices had two data blocks, one for all real components, one for all imaginary components -- "separate complex".)

I figure, since the storage now allows for it, that it should be possible to cast a complex-valued array to a real-valued array with twice as many elements, without copying the data.

MATLAB has a function typecast, which casts an array to a different type without copying data. It can be used, for example, to cast an array with 16 8-bit values to an array with 2 double floats. It does this without copying the data, the bit pattern is re-interpreted as the new type.

Sadly, this function does not work at all on complex-valued arrays.

I'm looking to replicate this code:

A = fftn(randn(40,60,20)); % some random complex-valued array
assert(~isreal(A))

sz = size(A);
B = reshape(A,1,[]);        % make into a vector
B = cat(1,real(B),imag(B)); % interleave real and imaginary values
B = reshape(B,[2,sz]);      % reshape back to original shape, with a new first dimension
assert(isreal(B))

The matrices A and B have (in R2018a and newer) the exact same data, in exactly the same order. However, to get to B we had to copy the data twice.

I tried creating a MEX-file that does this, but I don't see how to create a new array that references the data in the input array. This MEX-file works, but causes MATLAB to crash when clearing variables, because there are two arrays that reference the same data without them realizing that they share data (i.e. the reference count is not incremented).

// Build with:
//    mex -R2018a typecast_complextoreal.cpp

#include <mex.h>

#if MX_HAS_INTERLEAVED_COMPLEX==0
#error "This MEX-file must be compiled with the -R2018a flag"
#endif

#include <vector>

void mexFunction(int /*nlhs*/, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {

   // Validate input
   if(nrhs != 1) {
      mexErrMsgTxt("One input argument expected");
   }
   if(!mxIsDouble(prhs[0]) && !mxIsSingle(prhs[0])) {
    mexErrMsgTxt("Only floating-point arrays are supported");
   }

   // Get input array sizes
   mwSize nDims = mxGetNumberOfDimensions(prhs[0]);
   mwSize const* inSizes = mxGetDimensions(prhs[0]);

   // Create a 0x0 output matrix of the same type, but real-valued
   std::vector<mwSize> outSizes(nDims + 1, 0);
   plhs[0] = mxCreateNumericMatrix(0, 0, mxGetClassID(prhs[0]), mxREAL);

   // Set the output array data pointer to the input array's
   // NOTE! This is illegal, and causes MATLAB to crash when freeing both
   // input and output arrays, because it tries to free the same data
   // twice
   mxSetData(plhs[0], mxGetData(prhs[0]));

   // Set the output array sizes
   outSizes[0] = mxIsComplex(prhs[0]) ? 2 : 1;
   for(size_t ii = 0; ii < nDims; ++ii) {
      outSizes[ii + 1] = inSizes[ii];
   }
   mxSetDimensions(plhs[0], outSizes.data(), outSizes.size());
}

I'd love to hear of any ideas on how to proceed from here. I don't necessarily need to fix the MEX-file, if the solution is purely MATLAB code, so much the better.


Solution

  • See this FEX submission, which can do the complex --> 2 reals reinterpretation without copying the data (it can even point to interior contiguous sub-sections of the data without a copy):

    https://www.mathworks.com/matlabcentral/fileexchange/65842-sharedchild-creates-a-shared-data-copy-of-a-contiguous-subsection-of-an-existing-variable

    If you are simply reading and writing interleaved complex data files in R2018a and later, see this FEX submission:

    https://www.mathworks.com/matlabcentral/fileexchange/77530-freadcomplex-and-fwritecomplex