Higher-order Kronecker product
Since you are using Matlab, I suggest using SPGL1 as your reconstruction solver. It can handle complex numbers and you can specify your matrix operator as a function.

This example demonstrates how to specify your matrix operator as a function instead of a matrix (see opA).

In your case, I would define your operator as a function that:


takes a vector x as input,
reshapes that vector to a matrix X of your desired shape,
applies fft to the matrix: Z = fft(X) (Matlab's fft transforms each column separately if you give it a matrix),
applies fft to the transpose of the result: W = fft(Z'),
reshapes W back into a vector w and returns that


If your operation additionally includes some sub-sampling (which it usually does in order to be compressed sensing) then that can be implemented as indexing into w for the entries that are sampled.

The above steps imply that x is in the domain of the solution you are seeking and therefore sparse. Thus, the sparse vector here is in the time domain and you are sampling in the Fourier domain. If, on the other hand you want to model a signal that is sparse in frequency and you sample in the time domain, you should replace fft by ifft above.