MacResearch

Author: David Gohara

Continuing where we left off last time, we'll take our 2D FFT and convert it to a 2D planar FFT and then make a full-fledged 3D FFT out of it. This part of the tutorial should cruise by as we've covered most of the ground already.

You can download the Xcode project here

Let's start by setting our build target to part3 and selecting part3.c in the Groups and Files pane.

In order to do a pseudo-/full 3D FFT, we need to make a few changes. We've added a third dimension (nx), we are now taking the log2 of that dimension as well.

During the weights setup we expand that to include getting the MAX value across nx, ny and nz (although this isn't really necessary for part3):

fft_weights = create_fftsetup(MAX(MAX(logX,logY),logZ),kFFTRadix2);

3D FFT as 2D Planes

Some applications (such as EM image reconstruction) work on stacks of 2D images. This is similar to what we will do here. We define a new DSPSplitComplex type called tmp1 and will use that as a pointer to the correct location in memory of our data start for each plane.

int offset = 0;
for(i=0;i<nz;i++){
	tmp1.realp = input.realp + offset;
	tmp1.imagp = input.imagp + offset;
	
	fft2d_zip(fft_weights,&tmp1,1,0,logX,logY,FFT_FORWARD);
	
	offset += nx*ny;
}

Here we simply increment the pointers to the real and imaginary parts of the array by adding an offset value to pointer address at each loop iteration. Pretty simple. The rest of the code doesn't change.

Full 3D FFT

Continuing, set the build target to part4 and select part4.c

The next addition to the code is to define and allocate space for a results array. The results array will be used to rearrange the data so that we can perform the full 3D FFT. Recall that Accelerate doesn't include a 3D FFT by default (otherwise we would simply use that). So in order to perform one using Accelerate we need to first performa a 2D FFT, followed by a 1D FFT along the third dimension.

DSPSplitComplex input, result;
DSPSplitComplex tmp1;

input.realp = (float*)malloc(dims*sizeof(float));
input.imagp = (float*)malloc(dims*sizeof(float));

result.realp = (float*)malloc(dims*sizeof(float));
result.imagp = (float*)malloc(dims*sizeof(float));

...and of course at the end we include...

free(result.realp);
free(result.imagp);

Our 2D FFT will be performed on the planes as we did in part3. Once that is done we need to rearrange the data such that data along the Z-axis is now available in memory at a unit stride of 1. We can use another Accelerate function called mtrans (Matrix Transpose) to do this for us. mtrans is designed to work on 2D arrays, but we can achieve the same effect on a "3D" array by adjusting a few parameters:

mtrans(input.realp,1,result.realp,1,(nx*ny),nz);
mtrans(input.imagp,1,result.imagp,1,(nx*ny),nz);

input data array 	= input.realp (or input.imagp)
stride1 		= 1 (we want to process every element)
output data array 	= result.realp (the rearranged data)
stride2			= 1 (we want the data to be placed with a stride of 1)
dim1			= nx*ny (row count)
dim2			= nz (column count)

Most of this is straightforward. For dim1 we set that to nx*ny, because our row count is actually the full subset of one 2D plane. For dim2 our column count is the number of elements along the current Z direction. After the transpose these will be swapped (as you'll see below). Note that since our data is in split complex format we need to perform two of these operations (one for the real part and one for the imaginary part). mtrans (as it's engineered) will not work on interleaved complex data.

Once we have our rearranged data we can then perform a series of 1D FFT's. We will need to perform nx*ny FFT's each with a length of nz.

for(i=0;i<nx*ny;i++){
	tmp1.realp = &result.realp[i*nz];
	tmp1.imagp = &result.imagp[i*nz];
	
	fft_zip(fft_weights,&tmp1,1,logZ,FFT_FORWARD);
}

  We use fft_zip (1D complex in-place):

weights 	= fft_weights
data 		= tmp1 (which is a pointer to the mtrans results adjusted by i*nz)
signal stride 	= 1 (we want to process every element)
log2n 		= logZ (the base 2 exponent)
direction	= FFT_FORWARD (we are doing the forward transform)

Once all of the 1D FFT's have been completed, we need rearrange the data so that it's ordered the original way it was supplied. So we make another call (two actually) to mtrans. This time we will take the result data and put it back into the input data (overwriting what is there). We also swap the dim1 and dim2 values to ensure the data is placed properly:

mtrans(result.realp,1,input.realp,1,nz,(nx*ny));
mtrans(result.imagp,1,input.imagp,1,nz,(nx*ny));

We end up by converting to interleaved complex (ztoc) and then normalizing the data (vsmul) as before.

Checking the Result

To check that everything is internally consistent we can perform the forward and reverse FFTs to see if we get back our input values. Set the target to part5 and select part5.c. In both FFT loops we simply do a forward FFT followed by a reverse FFT.

Build and Run the project and you'll see in the run log that we get back what we put in.

To better understand what mtrans does and why it's such a powerhouse over simple copy operations (especially on Altivec), you may want to check out the following page on Apple's website:

http://developer.apple.com/hardwaredrivers/ve/algorithms.html#transpose