/*

    File:       vBigDSP.c

    Contains:   AltiVec-based Implementation of DSP routines (real & complex FFT, convolution)

    Version:    1.0

    Copyright:  © 1999 by Apple Computer, Inc., all rights reserved.

    Change History (most recent first):

                10/12/99    JK      Created

*/

/////////////////////////////////////////////////////////////////////////////////

//      File Name: vBigDSP.c                                                   //

//                                                                             //

//      This library provides a set of DSP routines, implemented using the     //

//      AltiVec instruction set.                                               //

//                                                                             //

//                                                                             //

//      Copyright © 1999 Apple Computer, Inc.  All rights reserved.            //

//                                                                             //

//      Version 1.0                                                            //

//                                                                             //                                                                             //

/////////////////////////////////////////////////////////////////////////////////

/*

    Disclaimer: IMPORTANT:  This Apple software is supplied to you by Apple Computer, Inc.

                ("Apple") in consideration of your agreement to the following terms, and your

                use, installation, modification or redistribution of this Apple software

                constitutes acceptance of these terms.  If you do not agree with these terms,

                please do not use, install, modify or redistribute this Apple software.

                In consideration of your agreement to abide by the following terms, and subject

                to these terms, Apple grants you a personal, non-exclusive license, under AppleÕs

                copyrights in this original Apple software (the "Apple Software"), to use,

                reproduce, modify and redistribute the Apple Software, with or without

                modifications, in source and/or binary forms; provided that if you redistribute

                the Apple Software in its entirety and without modifications, you must retain

                this notice and the following text and disclaimers in all such redistributions of

                the Apple Software.  Neither the name, trademarks, service marks or logos of

                Apple Computer, Inc. may be used to endorse or promote products derived from the

                Apple Software without specific prior written permission from Apple.  Except as

                expressly stated in this notice, no other rights or licenses, express or implied,

                are granted by Apple herein, including but not limited to any patent rights that

                may be infringed by your derivative works or by other works in which the Apple

                Software may be incorporated.

                The Apple Software is provided by Apple on an "AS IS" basis.  APPLE MAKES NO

                WARRANTIES, EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED

                WARRANTIES OF NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS FOR A PARTICULAR

                PURPOSE, REGARDING THE APPLE SOFTWARE OR ITS USE AND OPERATION ALONE OR IN

                COMBINATION WITH YOUR PRODUCTS.

                IN NO EVENT SHALL APPLE BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL OR

                CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE

                GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)

                ARISING IN ANY WAY OUT OF THE USE, REPRODUCTION, MODIFICATION AND/OR DISTRIBUTION

                OF THE APPLE SOFTWARE, HOWEVER CAUSED AND WHETHER UNDER THEORY OF CONTRACT, TORT

                (INCLUDING NEGLIGENCE), STRICT LIABILITY OR OTHERWISE, EVEN IF APPLE HAS BEEN

                ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

*/

#ifdef __MWERKS__

#include <AltiVec.h>

#endif

#include <errors.h>

#include <math.h>

#include <MacMemory.h>

#include "vBigDSP.h"

#ifndef PI

#define PI (3.1415926535897932384626433832795)

#endif

////////////////////////////////////////////////////////////////////////////////

//  temporary workspace buffer for FFT implementations 

////////////////////////////////////////////////////////////////////////////////

static Handle               gTempBufferHandle=nil;

////////////////////////////////////////////////////////////////////////////////

//  current buffer length for above handle

////////////////////////////////////////////////////////////////////////////////

static unsigned long        gTempBufferSize = 0;

////////////////////////////////////////////////////////////////////////////////

//  table of pre-calculated sin & cos values for use by ping pong FFT

////////////////////////////////////////////////////////////////////////////////

static Handle               gSinCosTableHandle = nil;

////////////////////////////////////////////////////////////////////////////////

//  length of currently allocated sin & cos table

////////////////////////////////////////////////////////////////////////////////

static unsigned long        gSinCosTableSize = 0;

////////////////////////////////////////////////////////////////////////////////

//

//  log2max

//

//  This function computes the ceiling of log2(n).  For example:

//

//  log2max(7) = 3

//  log2max(8) = 3

//  log2max(9) = 4

//

////////////////////////////////////////////////////////////////////////////////

static long log2max(long    n)

{

    long    power = 1;

    long    k = 1;

    if (n==1) {

        return 0;

    }

    while ((k <<= 1) < n) {

        power++;

    }

    return power;

}

////////////////////////////////////////////////////////////////////////////////

//

//  EnsureStaticBufferSize

//

// This function is used to ensure that the global Handle gTempBufferHandle

// is large enough to hold complexCount complex floats.  It compares the

// currently allocated length (if there is one) with the desired length,

// and, if necessary, allocates a new, larger handle.

// It will NOT shrink the handle, because a function higher in the calling

// chain may need the larger size.

////////////////////////////////////////////////////////////////////////////////

static OSErr EnsureStaticBufferSize(long    complexCount)

{

    OSErr       result = noErr;

    if (gTempBufferSize < complexCount) {

        gTempBufferSize = complexCount;

        if (gTempBufferHandle) {

            DisposeHandle(gTempBufferHandle);

            gTempBufferHandle = nil;

        }

        gTempBufferHandle = NewHandle(2*complexCount*sizeof(float));

        result = MemError();

        if (result != noErr) {

            gTempBufferHandle = nil;

            gTempBufferSize = 0;

        }

    }

    return result;

}

///////////////////////////////////////////////////////////////////////////////////////////////

// InitFFTSinCos

//  Initialize cos & sin lookup table.

///////////////////////////////////////////////////////////////////////////////////////////////

static OSErr InitFFTSinCos(long len) {

    OSErr           result;

    long            i;

    float           angle;

    float           *sinCosTable;

    float           curSin;

    float           curCos;

    long            handleSize;

    ////////////////////////////////////////////////////////////////////////////////

    // set length indicator to new length

    ////////////////////////////////////////////////////////////////////////////////

    gSinCosTableSize = len;

    ////////////////////////////////////////////////////////////////////////////////

    // deallocate old handle if there is one

    ////////////////////////////////////////////////////////////////////////////////

    if(gSinCosTableHandle) DisposeHandle(gSinCosTableHandle);

    ////////////////////////////////////////////////////////////////////////////////

    // calculate appropriate handle size for table.  Ensure a minimum length for

    //  a one-entry (two-value: cos & sin) table

    ////////////////////////////////////////////////////////////////////////////////

    handleSize = sizeof(float) * len;

    if (!handleSize)    handleSize = 2 * sizeof(float);

    ////////////////////////////////////////////////////////////////////////////////

    // alloc new handle for sin & cos table

    ////////////////////////////////////////////////////////////////////////////////

    gSinCosTableHandle = NewHandle(handleSize);

    result = MemError();

    if (result == noErr) {

        HLock(gSinCosTableHandle);

        ////////////////////////////////////////////////////////////////////////////////

        // get pointer to start of sin cos table

        // to reference as array

        ////////////////////////////////////////////////////////////////////////////////

        sinCosTable = *(float**)gSinCosTableHandle;

        if (len < 4) {

            ////////////////////////////////////////////////////////////////////////////////

            // always create first entries in table

            ////////////////////////////////////////////////////////////////////////////////            

            sinCosTable[0] = 1;

            sinCosTable[1] = 0;

            ////////////////////////////////////////////////////////////////////////////////

            // if len is 2, then add second pair to table

            ////////////////////////////////////////////////////////////////////////////////

            if (len == 2) {

                sinCosTable[2] = -1;

                sinCosTable[3] = 0;

            }

        } else {

            ////////////////////////////////////////////////////////////////////////////////

            // calculate cos & sin for range of 

            // [0, pi].

            ////////////////////////////////////////////////////////////////////////////////

            for(i=0; i < len/4; i++) {

                angle = 2.0 * PI * i/len;

                curCos = cos(angle);

                curSin = sin(angle);

                sinCosTable[2*i] = curCos;

                sinCosTable[2*i+1] = curSin;

                sinCosTable[len/2 + 2*i] = -curSin;

                sinCosTable[len/2 + 2*i + 1] = curCos;

            }

        }

        HUnlock(gSinCosTableHandle);

    }

    ////////////////////////////////////////////////////////////////////////////////

    // if there was an error, then dealloc

    // anything we allocated, and reset

    // length indicator to zero.

    ////////////////////////////////////////////////////////////////////////////////

    if (result != noErr) {

        if (gSinCosTableHandle != nil) {

            DisposeHandle(gSinCosTableHandle);

            gSinCosTableHandle = nil;

        }       

        gSinCosTableSize = 0;

    }

    return result;

}

///////////////////////////////////////////////////////////////////////////////////////////////

// ShutdownFFT

//

// deallocates any memory curently allocated for fft buffers.

///////////////////////////////////////////////////////////////////////////////////////////////

void ShutdownFFT()

{

    if (gSinCosTableHandle != nil) {

        DisposeHandle(gSinCosTableHandle);

        gSinCosTableHandle = nil;

    }

    gSinCosTableSize = 0;       

    if (gTempBufferHandle != nil) {

        DisposeHandle(gTempBufferHandle);

        gTempBufferHandle = nil;    

    }

    gTempBufferSize = 0;

}

#pragma mark -

#pragma mark R E C U R S I V E   C O M P L E X

////////////////////////////////////////////////////////////////////////////////

//  SquareComplexTransposeVector performs a transpose on a square matrix of

//  complex floats with side length of rowLength.  Data is assumed to be

//  arranged lexicographically, i.e., with side length n:

//

//  re[0] im[0] re[1] im[1] ... re[n-1] im[n-1]

//  re[n] im[n]     ...         re[2n-1] im[2n-1]

//  .

//  .

//  re[(n-1)*n] im[(n-1)*n] ... re[n^2-1] im[n^2-1]

//  

//  Data becomes:

//

//  re[0] im[0] re[n] im[n] re[2n] im[2n]...

//  re[1] im[1] re[n+1] im[n+1] ...

//  re[2] im[2] ...

//  .

//  .

//  .

//  re[n-1] im[n-1] ...     re[n^2-1] im[n^2-1]

//

//

////////////////////////////////////////////////////////////////////////////////

static void SquareComplexTransposeVector(float *data, long rowLength)

{

    long                    i, j;

    vector float            *pInLeft1, *pInLeft2;

    vector float            *pInTop1, *pInTop2;

    vector float            vInLeft1, vInLeft2;

    vector float            vInTop1, vInTop2;

    vector float            vOutLeft1, vOutLeft2;

    vector float            vOutTop1, vOutTop2;

    vector unsigned char    vMergeHiPairPerm;

    vector unsigned char    vMergeLoPairPerm;       

    ////////////////////////////////////////////////////////////////////////////////

    // initialize a permute vector that, given input vectors:

    //

    //  X = x0 x1 x2 x3

    //  Y = y0 y1 y2 y3

    //

    // will generate

    //  Z = x0 x1 y0 y1

    ////////////////////////////////////////////////////////////////////////////////

    vMergeHiPairPerm = (vector unsigned char)(0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23);

    ////////////////////////////////////////////////////////////////////////////////

    // initialize a permute vector that, given input vectors:

    //

    //  X = x0 x1 x2 x3

    //  Y = y0 y1 y2 y3

    //

    // will generate

    //  Z = x2 x3 y2 y3

    ////////////////////////////////////////////////////////////////////////////////

    vMergeLoPairPerm = (vector unsigned char)(8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31);

    for (i=0; i<(rowLength/2); i++) {

        ////////////////////////////////////////////////////////////////////////////////

        // set pointers to point to the beginning of row 2i and 2i+1

        ////////////////////////////////////////////////////////////////////////////////

        pInLeft1 = ((vector float*)data) + i * rowLength;

        pInLeft2 = pInLeft1 + rowLength / 2;

        ////////////////////////////////////////////////////////////////////////////////

        // set pointers to point to the first elements of columns 2i and 2i+1

        ////////////////////////////////////////////////////////////////////////////////

        pInTop1 = ((vector float*)data) + i;

        pInTop2 = pInTop1 + rowLength / 2;

        for (j=0; j<i; j++) {

            ////////////////////////////////////////////////////////////////////////////////

            // read in two vectors that contain a 2x2 square of matrix elements

            // from the left side of the diagonal

            ////////////////////////////////////////////////////////////////////////////////

            vInLeft1 = *pInLeft1;

            vInLeft2 = *pInLeft2;

            ////////////////////////////////////////////////////////////////////////////////

            // read in two vectors that contain a 2x2 square of matrix elements

            // from the top side of the diagonal

            ////////////////////////////////////////////////////////////////////////////////

            vInTop1 = *pInTop1;

            vInTop2 = *pInTop2;

            ////////////////////////////////////////////////////////////////////////////////

            //  We now transpose our two 2x2 sub-matrices, and swap their positions.  The

            // transpose is done with vector permute operations.

            ////////////////////////////////////////////////////////////////////////////////

            vOutLeft1 = vec_perm(vInTop1, vInTop2, vMergeHiPairPerm);

            vOutLeft2 = vec_perm(vInTop1, vInTop2, vMergeLoPairPerm);

            vOutTop1 = vec_perm(vInLeft1, vInLeft2, vMergeHiPairPerm);

            vOutTop2 = vec_perm(vInLeft1, vInLeft2, vMergeLoPairPerm);

            ////////////////////////////////////////////////////////////////////////////////

            // store the transposed matrices in their swapped positions

            ////////////////////////////////////////////////////////////////////////////////

            *pInLeft1 = vOutLeft1;

            *pInLeft2 = vOutLeft2;

            *pInTop1 = vOutTop1;

            *pInTop2 = vOutTop2;

            ////////////////////////////////////////////////////////////////////////////////

            // advance pointers to next elements in the current columns & rows

            ////////////////////////////////////////////////////////////////////////////////

            pInLeft1 += 1;

            pInLeft2 += 1;

            pInTop1 += rowLength;

            pInTop2 += rowLength;

        }

        ////////////////////////////////////////////////////////////////////////////////

        // when the above loop is finished, our input pointers point to the 2x2 sub-

        // matrix that is on the diagonal of our matrix to be transposed.  We read in

        // this sub-matrix into two vectors, transpose the sub-matrix (using vector

        // transposes) and store the sub-matrix.

        ////////////////////////////////////////////////////////////////////////////////

        vInLeft1 = *pInLeft1;

        vInLeft2 = *pInLeft2;

        vOutLeft1 = vec_perm(vInLeft1, vInLeft2, vMergeHiPairPerm);

        vOutLeft2 = vec_perm(vInLeft1, vInLeft2, vMergeLoPairPerm);

        *pInLeft1 = vOutLeft1;

        *pInLeft2 = vOutLeft2;

    }

}

////////////////////////////////////////////////////////////////////////////////

//  SquareComplexTransposeTwist performs a transpose on a square matrix of

//  complex floats with side length of rowLength.  Data is assumed to be

//  arranged lexicographically, i.e., with side length n.

//

//

//

//  re[0] im[0] re[1] im[1] ... re[n-1] im[n-1]

//  re[n] im[n]     ...         re[2n-1] im[2n-1]

//  .

//  .

//  re[(n-1)*n] im[(n-1)*n] ... re[n^2-1] im[n^2-1]

//  

//

//  In addition to performing the transpose, each element X(j, k) is

//  multiplied by the twist factor e^(+/- 2 pi i j k / (n^2) ). 

//

////////////////////////////////////////////////////////////////////////////////

static void SquareComplexTransposeTwist(float *data, long sideLength, long isInverse)

{

    long                        i, j;

    vector float            *pInLeft1, *pInLeft2;

    vector float            *pInTop1, *pInTop2;

    vector float            vInLeft1, vInLeft2;

    vector float            vInTop1, vInTop2;

    vector float            vOutLeft1, vOutLeft2;

    vector float            vOutTop1, vOutTop2;

    vector float            vTwistedTop1, vTwistedTop2;

    vector float            vPermutedLeft1, vPermutedLeft2;

    vector unsigned char    vMergeHiPairPerm;

    vector unsigned char    vMergeLoPairPerm;       

    vector unsigned char    vSwappedPerm;

    vector float            vCosTwistA0;

    vector float            vCosTwistA1;

    vector float            vCosTemp0, vCosTemp1;

    vector float            vSinTwistA0;

    vector float            vSinTwistA1;

    vector float            vA;

    vector float            vB;

    vector float            vTransition;

    vector float            vSinSignMultiplier;

    vector float            vZero = (vector float)(0);

    double                  fTemp;

    double                  baseAngle1;

    double                  baseAngle2;

    double                  fSinBaseAngle1;

    double                  fSinBaseAngle2;

    ////////////////////////////////////////////////////////////////////////////////

    // initialize multiplier for sin vector, for whether we are multiplying

    // by e^(2 pi i l m / (n^2) ) or e^(-2 pi i l m / (n^2) ).

    ////////////////////////////////////////////////////////////////////////////////

    if (isInverse) {

        vSinSignMultiplier = (vector float)(-1, 1, -1, 1);

    } else {

        vSinSignMultiplier = (vector float)(1, -1, 1, -1);  

    }

    ////////////////////////////////////////////////////////////////////////////////

    // initialize a permute vector that, given input vectors:

    //

    //  X = x0 x1 x2 x3

    //  Y = y0 y1 y2 y3

    //

    // will generate

    //  Z = x0 x1 y0 y1

    ////////////////////////////////////////////////////////////////////////////////

    vMergeHiPairPerm = (vector unsigned char)(0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23);

    ////////////////////////////////////////////////////////////////////////////////

    // initialize a permute vector that, given input vectors:

    //

    //  X = x0 x1 x2 x3

    //  Y = y0 y1 y2 y3

    //

    // will generate

    //  Z = x2 x3 y2 y3

    ////////////////////////////////////////////////////////////////////////////////

    vMergeLoPairPerm = (vector unsigned char)(8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31);

    ////////////////////////////////////////////////////////////////////////////////

    // initialize a permute vector that, given input vectors:

    //

    //  X = x0 x1 x2 x3

    //  Y = y0 y1 y2 y3

    //

    // will generate

    //  Z = x1 x0 x3 y2

    ////////////////////////////////////////////////////////////////////////////////

    vSwappedPerm  = (vector unsigned char)(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11);

    ////////////////////////////////////////////////////////////////////////////////

    //

    // special-case upper-left 2x2 sub-matrix

    //

    ////////////////////////////////////////////////////////////////////////////////

    ////////////////////////////////////////////////////////////////////////////////

    //  set pointers to point at first and second rows  

    ////////////////////////////////////////////////////////////////////////////////

    pInLeft1 = ((vector float*)data);

    pInLeft2 = pInLeft1 + sideLength / 2;

    ////////////////////////////////////////////////////////////////////////////////

    //  create twist multiplier cos & sin vectors

    //

    //  the vTransition vector is used to move elements from the scalar domain

    //  to the vector domain.

    ////////////////////////////////////////////////////////////////////////////////

    baseAngle1 = 0;

    baseAngle2 = (2*PI)/(sideLength*sideLength);

    fSinBaseAngle1 = sin(baseAngle1);

    fSinBaseAngle2 = sin(baseAngle2);

    ((float*)&vTransition)[0] = 1;

    ((float*)&vTransition)[2] = cos(baseAngle1);

    ((float*)&vTransition)[3] = cos(baseAngle2);

    vCosTwistA0 = vec_splat(vTransition, 0);

    vCosTwistA1 = vec_mergel(vTransition, vTransition);

    vSinTwistA0 = vZero;

    ((float*)&vTransition)[2] = fSinBaseAngle1;

    ((float*)&vTransition)[3] = fSinBaseAngle2;

    vSinTwistA1 = vec_mergel(vTransition, vTransition);

    vSinTwistA1 = vec_madd(vSinTwistA1, vSinSignMultiplier, vZero);

    ////////////////////////////////////////////////////////////////////////////////

    // twist multiply & permute 2x2 square on diagonal.

    ////////////////////////////////////////////////////////////////////////////////

    vInLeft1 = *pInLeft1;

    vInLeft2 = *pInLeft2;

    vPermutedLeft1 = vec_perm(vInLeft1, vInLeft2, vMergeHiPairPerm);

    vPermutedLeft2 = vec_perm(vInLeft1, vInLeft2, vMergeLoPairPerm);

    vOutLeft1 = vec_madd(vPermutedLeft1, vCosTwistA0, vZero);

    vOutLeft1 = vec_madd(vec_perm(vPermutedLeft1, vPermutedLeft1, vSwappedPerm), vSinTwistA0, vOutLeft1);

    vOutLeft2 = vec_madd(vPermutedLeft2, vCosTwistA1, vZero);

    vOutLeft2 = vec_madd(vec_perm(vPermutedLeft2, vPermutedLeft2, vSwappedPerm), vSinTwistA1, vOutLeft2);

    *pInLeft1 = vOutLeft1;

    *pInLeft2 = vOutLeft2;

    ////////////////////////////////////////////////////////////////////////////////

    // loop through all remaining rows/columns 2 at a time

    ////////////////////////////////////////////////////////////////////////////////

    for (i=1; i<(sideLength/2); i++) {

        ////////////////////////////////////////////////////////////////////////////////

        // set pointers to point to the beginning of row 2i and 2i+1

        ////////////////////////////////////////////////////////////////////////////////

        pInLeft1 = ((vector float*)data) + i * sideLength;

        pInLeft2 = pInLeft1 + sideLength / 2;

        ////////////////////////////////////////////////////////////////////////////////

        // set pointers to point to the first elements of columns 2i and 2i+1

        ////////////////////////////////////////////////////////////////////////////////

        pInTop1 = ((vector float*)data) + i;

        pInTop2 = pInTop1 + sideLength / 2;

        ////////////////////////////////////////////////////////////////////////////////

        // set up vectors for incremental updating of cos & sin vectors

        //

        // Given c = cos(w) and s = sin(w), if we want to find the cos and sin of w plus

        // some small angle d, then we can do so by defining:

        //  

        // a = 2 * ((sin(d/2)) ^ 2)

        // b = sin(d)

        // 

        // Then, we can calculate the cos and sin of the updated angle w+d as:

        // 

        // cos(w+d) = c - ac - bs

        // sin(w+d) = s - as + bc

        // 

        // the vectors vA and vB contain these increment factors a and b.  Note that

        // each outer loop iteration requires different increment angles.  Looking at

        // the matrix, the angle multipliers that we use for our cos & sin twist

        // multipliers look something like:

        //

        //      0/N     0/N     0/N     0/N     0/N     0/N     0/N ... 0/N

        //      0/N     1/N     2/N     3/N     4/N     5/N     6/N     

        //      0/N     2/N     4/N     6/N     8/N     10/N    12/N    .

        //      0/N     3/N     6/N     9/N     12/N    15/N    18/N    .   

        //      0/N     4/N     8/N     12/N    16/N    20/N    24/N

        //      .       .

        //      .       .

        //      0/N                                         ...         (N-1)(N-1)

        //

        // Since we're doing two rows/columns at a time, the vA and vB vectors

        // need to have two different sets of increment values.  For example, for

        // the third and fourth columns (handled in the second time through the

        // outer loop) each row increments by 2/N and 3/N.  Since we actually 

        // have two sets of twist multiplier vectors, we skip a row for each, so

        // the increment angle is doubled.  This helps reduce calculation dependencies

        // and cuts the number of increment iterations in half, which helps reduce

        // errors in the single-precision calculations used.

        //      

        ////////////////////////////////////////////////////////////////////////////////

        baseAngle1 = 4*i*PI/(sideLength*sideLength);

        baseAngle2 = (2*(2*i+1)*PI)/(sideLength*sideLength);

        fSinBaseAngle1 = sin(baseAngle1);

        fTemp = 2*fSinBaseAngle1*fSinBaseAngle1;

        ((float*)&vTransition)[0] = fTemp;

        ((float*)&vTransition)[2] = sin(2*baseAngle1);

        fSinBaseAngle2 = sin(baseAngle2);

        fTemp = 2*fSinBaseAngle2*fSinBaseAngle2;

        ((float*)&vTransition)[1] = fTemp;

        ((float*)&vTransition)[3] = sin(2*baseAngle2);

        vA = vec_mergeh(vTransition, vTransition);

        vB = vec_mergel(vTransition, vTransition);

        vB = vec_madd(vB, vSinSignMultiplier, vZero);

        ////////////////////////////////////////////////////////////////////////////////

        //  Set up the initial twist multiplier vectors for the beginning of this

        // row/column pair.  Values are calculated in the scalar domain, and then

        // transferred to the vector domain via the vTransition vector.

        ////////////////////////////////////////////////////////////////////////////////

        ((float*)&vTransition)[0] = 1;

        ((float*)&vTransition)[2] = cos(baseAngle1);

        ((float*)&vTransition)[3] = cos(baseAngle2);

        vCosTwistA0 = vec_splat(vTransition, 0);

        vCosTwistA1 = vec_mergel(vTransition, vTransition);

        vSinTwistA0 = vZero;

        ((float*)&vTransition)[2] = fSinBaseAngle1;

        ((float*)&vTransition)[3] = fSinBaseAngle2;

        vSinTwistA1 = vec_mergel(vTransition, vTransition);

        vSinTwistA1 = vec_madd(vSinTwistA1, vSinSignMultiplier, vZero);

        ////////////////////////////////////////////////////////////////////////////////

        // load in the leftmost and topmost 2x2 matrices for our row/column pair

        ////////////////////////////////////////////////////////////////////////////////

        vInLeft1 = *pInLeft1;

        vInLeft2 = *pInLeft2;

        vInTop1 = *pInTop1;

        vInTop2 = *pInTop2;

        ////////////////////////////////////////////////////////////////////////////////

        // twist multiply in top vectors and use permute to transpose before storing

        // to position in new row.

        ////////////////////////////////////////////////////////////////////////////////

        vTwistedTop1 = vInTop1;

        vTwistedTop2 = vec_madd(vInTop2, vCosTwistA1, vZero);

        vTwistedTop2 = vec_madd(vec_perm(vInTop2, vInTop2, vSwappedPerm), vSinTwistA1, vTwistedTop2);

        vOutLeft1 = vec_perm(vTwistedTop1, vTwistedTop2, vMergeHiPairPerm);

        vOutLeft2 = vec_perm(vTwistedTop1, vTwistedTop2, vMergeLoPairPerm);

        ////////////////////////////////////////////////////////////////////////////////

        // use permute to transpose 2x2 matrix from current row, and then twist

        // multiply.

        ////////////////////////////////////////////////////////////////////////////////

        vPermutedLeft1 = vec_perm(vInLeft1, vInLeft2, vMergeHiPairPerm);

        vPermutedLeft2 = vec_perm(vInLeft1, vInLeft2, vMergeLoPairPerm);

        vOutTop1 = vPermutedLeft1;

        vOutTop2 = vec_madd(vPermutedLeft2, vCosTwistA1, vZero);

        vOutTop2 = vec_madd(vec_perm(vPermutedLeft2, vPermutedLeft2, vSwappedPerm), vSinTwistA1, vOutTop2);

        ////////////////////////////////////////////////////////////////////////////////

        // update sin & cos twist multiplier vectors

        ////////////////////////////////////////////////////////////////////////////////

        vCosTemp0 = vec_nmsub(vCosTwistA0, vA, vCosTwistA0);

        vCosTemp0 = vec_nmsub(vSinTwistA0, vB, vCosTemp0);

        vSinTwistA0 = vec_nmsub(vSinTwistA0, vA, vSinTwistA0);

        vSinTwistA0 = vec_madd (vCosTwistA0, vB, vSinTwistA0);

        vCosTwistA0 = vCosTemp0;

        vCosTemp1 = vec_nmsub(vCosTwistA1, vA, vCosTwistA1);

        vCosTemp1 = vec_nmsub(vSinTwistA1, vB, vCosTemp1);

        vSinTwistA1 = vec_nmsub(vSinTwistA1, vA, vSinTwistA1);

        vSinTwistA1 = vec_madd (vCosTwistA1, vB, vSinTwistA1);

        vCosTwistA1 = vCosTemp1;

        ////////////////////////////////////////////////////////////////////////////////

        // store twist multiplied and transposed matrices.

        ////////////////////////////////////////////////////////////////////////////////

        *pInLeft1 = vOutLeft1;

        *pInLeft2 = vOutLeft2;

        *pInTop1 = vOutTop1;

        *pInTop2 = vOutTop2;

        ////////////////////////////////////////////////////////////////////////////////

        // update pointers to move down in current columns and right in current rows

        ////////////////////////////////////////////////////////////////////////////////

        pInLeft1 += 1;

        pInLeft2 += 1;

        pInTop1 += sideLength;

        pInTop2 += sideLength;

        ////////////////////////////////////////////////////////////////////////////////

        // loop through remaining 2x2 sub-matrices in current row/column pair until

        // we reach the diagonal of the matrix.     

        ////////////////////////////////////////////////////////////////////////////////

        for (j=1; j<i; j++) {

            ////////////////////////////////////////////////////////////////////////////////

            // load in 2x2 matrices for our row/column pair

            ////////////////////////////////////////////////////////////////////////////////

            vInLeft1 = *pInLeft1;

            vInLeft2 = *pInLeft2;

            vInTop1 = *pInTop1;

            vInTop2 = *pInTop2;

            ////////////////////////////////////////////////////////////////////////////////

            // twist multiply in top vectors and use permute to transpose before storing

            // to position in new row.

            ////////////////////////////////////////////////////////////////////////////////

            vTwistedTop1 = vec_madd(vInTop1, vCosTwistA0, vZero);

            vTwistedTop1 = vec_madd(vec_perm(vInTop1, vInTop1, vSwappedPerm), vSinTwistA0, vTwistedTop1);

            vTwistedTop2 = vec_madd(vInTop2, vCosTwistA1, vZero);

            vTwistedTop2 = vec_madd(vec_perm(vInTop2, vInTop2, vSwappedPerm), vSinTwistA1, vTwistedTop2);

            vOutLeft1 = vec_perm(vTwistedTop1, vTwistedTop2, vMergeHiPairPerm);

            vOutLeft2 = vec_perm(vTwistedTop1, vTwistedTop2, vMergeLoPairPerm);

            ////////////////////////////////////////////////////////////////////////////////

            // use permute to transpose 2x2 matrix from current row, and then twist

            // multiply.

            ////////////////////////////////////////////////////////////////////////////////

            vPermutedLeft1 = vec_perm(vInLeft1, vInLeft2, vMergeHiPairPerm);

            vPermutedLeft2 = vec_perm(vInLeft1, vInLeft2, vMergeLoPairPerm);

            vOutTop1 = vec_madd(vPermutedLeft1, vCosTwistA0, vZero);

            vOutTop1 = vec_madd(vec_perm(vPermutedLeft1, vPermutedLeft1, vSwappedPerm), vSinTwistA0, vOutTop1);

            vOutTop2 = vec_madd(vPermutedLeft2, vCosTwistA1, vZero);

            vOutTop2 = vec_madd(vec_perm(vPermutedLeft2, vPermutedLeft2, vSwappedPerm), vSinTwistA1, vOutTop2);

            ////////////////////////////////////////////////////////////////////////////////

            // update sin & cos twist multiplier vectors

            ////////////////////////////////////////////////////////////////////////////////

            vCosTemp0 = vec_nmsub(vCosTwistA0, vA, vCosTwistA0);

            vCosTemp0 = vec_nmsub(vSinTwistA0, vB, vCosTemp0);

            vSinTwistA0 = vec_nmsub(vSinTwistA0, vA, vSinTwistA0);

            vSinTwistA0 = vec_madd (vCosTwistA0, vB, vSinTwistA0);

            vCosTwistA0 = vCosTemp0;

            vCosTemp1 = vec_nmsub(vCosTwistA1, vA, vCosTwistA1);

            vCosTemp1 = vec_nmsub(vSinTwistA1, vB, vCosTemp1);

            vSinTwistA1 = vec_nmsub(vSinTwistA1, vA, vSinTwistA1);

            vSinTwistA1 = vec_madd (vCosTwistA1, vB, vSinTwistA1);

            vCosTwistA1 = vCosTemp1;

            ////////////////////////////////////////////////////////////////////////////////

            // store twist multiplied and transposed matrices.

            ////////////////////////////////////////////////////////////////////////////////

            *pInLeft1 = vOutLeft1;

            *pInLeft2 = vOutLeft2;

            *pInTop1 = vOutTop1;

            *pInTop2 = vOutTop2;

            ////////////////////////////////////////////////////////////////////////////////

            // update pointers to move down in current columns and right in current rows

            ////////////////////////////////////////////////////////////////////////////////

            pInLeft1 += 1;

            pInLeft2 += 1;

            pInTop1 += sideLength;

            pInTop2 += sideLength;

        }

        ////////////////////////////////////////////////////////////////////////////////

        // At this point, the above loop has left the pointers at the 2x2 sub-matrix

        // on the diagonal.  We read in this 2x2 matrix, transpose it, twist multiply

        // it, and store it back.

        ////////////////////////////////////////////////////////////////////////////////

        vInLeft1 = *pInLeft1;

        vInLeft2 = *pInLeft2;

        vPermutedLeft1 = vec_perm(vInLeft1, vInLeft2, vMergeHiPairPerm);

        vPermutedLeft2 = vec_perm(vInLeft1, vInLeft2, vMergeLoPairPerm);

        vOutLeft1 = vec_madd(vPermutedLeft1, vCosTwistA0, vZero);

        vOutLeft1 = vec_madd(vec_perm(vPermutedLeft1, vPermutedLeft1, vSwappedPerm), vSinTwistA0, vOutLeft1);

        vOutLeft2 = vec_madd(vPermutedLeft2, vCosTwistA1, vZero);

        vOutLeft2 = vec_madd(vec_perm(vPermutedLeft2, vPermutedLeft2, vSwappedPerm), vSinTwistA1, vOutLeft2);

        *pInLeft1 = vOutLeft1;

        *pInLeft2 = vOutLeft2;

    }

}

////////////////////////////////////////////////////////////////////////////////

//  fft_recursive

//

//  Performs a recursive (N/2 by 2) forward or inverse FFT on the

// signal pointed to by pData. This is done in the following manner.  First, the

// signal X is divided into two parts, namely X1 = X(0, N/2-1) and 

// X2 = X(N/2, N-1). These two sub-signals are then turned into the sum and

// difference, respectively, of these two signals. This is in effect an FFT

// of length 2 performed between each pair of elements taken from X1 and X2.

// This produces two new subsignals X1' and X2'.

//

// Next, a twist operation is performed on X2', where X2'(j) is multiplied

// by e^(2 pi ij / N).

//

// Next, a length N/2 complex FFT is performed on each of X1' and X2', yielding

// X1" and X2".

//

// Finally, the elements of X1" and X2" are merged, so that the resulting

// length-N signal Y is given by:

//

//  Y = X1"(0) X2"(0) X1"(1) X2"(1) ... X1"(N/2-1) X2"(N/2-1)

//

//  This resulting signal Y is the FFT of the original signal X

//

////////////////////////////////////////////////////////////////////////////////

static OSErr fft_recursive(float *pData, unsigned long len, long isign)

{

    OSErr                       result;

    long                        i;

    vector float                *pInLoBottom, *pInLoTop, *pInHiBottom, *pInHiTop;

    vector float                *pOutLoBottom, *pOutLoTop, *pOutHiBottom, *pOutHiTop;

    vector float                vLoBottom, vLoTop, vHiBottom, vHiTop;

    vector float                vZero;

    vector float                vSinLo;

    vector float                vCosLo;

    vector float                vSinHi;

    vector float                vCosHi;

    vector float                vFTransition;

    vector float                vFNegTransition;

    double                      updateA;

    double                      updateB;

    double                      newCos1;

    double                      newCos2;

    double                      newSin1;

    double                      newSin2;

    double                      tempCos1;

    double                      tempCos2;

    double                      startSinMul;

    vector unsigned char        vCosLoPerm;

    vector unsigned char        vSinPerm;

    vector unsigned char        vSwapPairPerm;

    vector unsigned char        vMergeLoPairPerm;

    vector unsigned char        vMergeHiPairPerm;

    vector unsigned long        vHiPairSelect;

    vector float                vDiffBottom;

    vector float                vDiffTop;   

    vector float                vOneHalf;

    vector float                vTwistBottom;

    vector float                vTwistTop;

    vector float                vSwappedDiffBottom;

    vector float                vSwappedDiffTop;

    ////////////////////////////////////////////////////////////////////////////////

    // make sure that work buffer is large enough to hold half of entire data

    // length for merge operation at end of routine.

    ////////////////////////////////////////////////////////////////////////////////

    result = EnsureStaticBufferSize(len/2);

    if (result == noErr) {

        ////////////////////////////////////////////////////////////////////////////////

        // initialize a permute vector that, given input vectors:

        //

        //  X = x0 x1 x2 x3

        //  Y = y0 y1 y2 y3

        //

        // will generate

        //  Z = x1 x0 x3 x2

        ////////////////////////////////////////////////////////////////////////////////

        vSwapPairPerm   = (vector unsigned char)(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11);

        ////////////////////////////////////////////////////////////////////////////////

        // initialize a permute vector that, given input vectors:

        //

        //  X = x0 x1 x2 x3

        //  Y = y0 y1 y2 y3

        //

        // will generate

        //  Z = x0 x0 y3 y3

        ////////////////////////////////////////////////////////////////////////////////

        vCosLoPerm      = (vector unsigned char)(0, 1, 2, 3, 0, 1, 2, 3, 28, 29, 30, 31, 28, 29, 30, 31);

        ////////////////////////////////////////////////////////////////////////////////

        // initialize a select vector that, given input vectors:

        //

        //  X = x0 x1 x2 x3

        //  Y = y0 y1 y2 y3

        //

        // will generate

        //  Z = x0 x1 y2 y3

        ////////////////////////////////////////////////////////////////////////////////

        vHiPairSelect   = (vector unsigned long)((vector signed long)(0, 0, -1, -1));

        if (isign == -1) {

            startSinMul = 1;

            ////////////////////////////////////////////////////////////////////////////////

            // initialize a permute vector that, given input vectors:

            //

            //  X = x0 x1 x2 x3

            //  Y = y0 y1 y2 y3

            //

            // will generate

            //  Z = x0 y0 x1 y1

            ////////////////////////////////////////////////////////////////////////////////

            vSinPerm = (vector unsigned char)(0, 1, 2, 3, 16, 17, 18, 19, 4, 5, 6, 7, 20, 21, 22, 23);

        } else {

            startSinMul = -1;

            ////////////////////////////////////////////////////////////////////////////////

            // initialize a permute vector that, given input vectors:

            //

            //  X = x0 x1 x2 x3

            //  Y = y0 y1 y2 y3

            //

            // will generate

            //  Z = y0 x0 y1 x1

            ////////////////////////////////////////////////////////////////////////////////

            vSinPerm = (vector unsigned char)(16, 17, 18, 19, 0, 1, 2, 3, 20, 21, 22, 23, 4, 5, 6, 7);      

        }

        ////////////////////////////////////////////////////////////////////////////////

        // initialize a permute vector that, given input vectors:

        //

        //  X = x0 x1 x2 x3

        //  Y = y0 y1 y2 y3

        //

        // will generate

        //  Z = x0 x1 y0 y1

        ////////////////////////////////////////////////////////////////////////////////

        vMergeHiPairPerm = (vector unsigned char)(0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23);

        ////////////////////////////////////////////////////////////////////////////////

        // initialize a permute vector that, given input vectors:

        //

        //  X = x0 x1 x2 x3

        //  Y = y0 y1 y2 y3

        //

        // will generate

        //  Z = x2 x3 y2 y3

        ////////////////////////////////////////////////////////////////////////////////

        vMergeLoPairPerm = (vector unsigned char)(8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31);

        vZero = (vector float)(0);

        ////////////////////////////////////////////////////////////////////////////////

        // start source pointers at start and end of high and low halves of the data

        ////////////////////////////////////////////////////////////////////////////////

        pInLoBottom = (vector float*)pData;

        pInHiBottom = pInLoBottom + len/4;

        pInLoTop = pInHiBottom - 1;

        pInHiTop = pInLoBottom + (len/2) - 1;

        ////////////////////////////////////////////////////////////////////////////////

        // start dest pointers at start and end of high and low halves of the data

        ////////////////////////////////////////////////////////////////////////////////

        pOutHiBottom = pInHiBottom;

        pOutHiTop = pInHiTop;

        pOutLoBottom = pInLoBottom;

        pOutLoTop = (pOutLoBottom + len/4) - 1;

        ////////////////////////////////////////////////////////////////////////////////

        //

        //  Given c = cos(w) and s = sin(w), if we

        //  want to find the cos and sin of w plus

        //  some small angle d, then we can do so

        //  by defining:

        //  

        //  a = 2 * ((sin(d/2)) ^ 2)

        //  b = sin(d)

        //  

        //  Then, we can calculate the cos and sin

        //  of the updated angle w+d as:

        //  

        //  cos(w+d) = c - ac - bs

        //  sin(w+d) = s - as + bc

        //  

        //  We will need to incrementally calculate

        //  cos(w) and sin(w), so we define the

        //  following scalar values.  We use scalar

        //  because we need the precision of doubles,

        //  and vectors are floats. 

        //

        ////////////////////////////////////////////////////////////////////////////////

        updateA = sin(2*PI/len);

        updateA *= updateA*2;

        updateB = sin(4*PI/len);

        ////////////////////////////////////////////////////////////////////////////////

        //

        // We want to have four cos and sin vectors

        // that are updated incrementally, using

        // the doubles that we have used to calculate

        // our new values.  To do this we have a

        // "transition" vector that allows us to get

        // our scalar values into the vector domain.

        //

        // Scalar values are stored into the elements

        // of the transition vector, and then 

        // our end cos and sin vectors are created

        // by permuting values out of our transition

        // vector and other previously created sin

        // and cos vectors.

        // 

        ////////////////////////////////////////////////////////////////////////////////

        ////////////////////////////////////////////////////////////////////////////////

        // set up vSinLo as

        // vSinLo = -sin(0) sin(0) -sin(pi/len) sin(pi/len)

        ////////////////////////////////////////////////////////////////////////////////

        ((float*)&vFTransition)[0] = 0;

        ((float*)&vFTransition)[1] = 0;

        newSin1 = sin(2*PI/len);

        ((float*)&vFTransition)[2] = newSin1*startSinMul;

        ((float*)&vFTransition)[3] = -newSin1*startSinMul;

        vSinLo = vFTransition;

        ////////////////////////////////////////////////////////////////////////////////

        // set up vSinHi as

        // vSinHi = -sin((len-2)*pi/len) sin((len-2)*pi/len) -sin((len-1)*pi/len) sin((len-1)*pi/len)

        ////////////////////////////////////////////////////////////////////////////////

        newSin2 = sin(4*PI/len);

        ((float*)&vFTransition)[0] = newSin2*startSinMul;   

        ((float*)&vFTransition)[1] = -newSin2*startSinMul;

        vSinHi = vFTransition;

        ///////////////////////////////////////////// 

        // set up vCosLo as

        // vCosLo = -cos(0) cos(0) -cos(pi/len) cos(pi/len)

        ////////////////////////////////////////////////////////////////////////////////

        ((float*)&vFTransition)[0] = 1;

        ((float*)&vFTransition)[1] = 1;

        newCos1 = cos(2*PI/len);

        ((float*)&vFTransition)[2] = newCos1;

        ((float*)&vFTransition)[3] = newCos1;

        vCosLo = vFTransition;

        ////////////////////////////////////////////////////////////////////////////////

        // set up vCosHi as

        // vCosHi = -cos((len-2)*pi/len) cos((len-2)*pi/len) -cos((len-1)*pi/len) cos((len-1)*pi/len)

        ////////////////////////////////////////////////////////////////////////////////

        newCos2 = cos(4*PI/len);

        ((float*)&vFTransition)[0] = newCos2;   

        ((float*)&vFTransition)[1] = newCos2;

        vCosHi = vec_sub(vZero, vFTransition);  

        if (isign == -1) {

            ////////////////////////////////////////////////////////////////////////////////

            // we are doing a forward FFT

            ////////////////////////////////////////////////////////////////////////////////

            for (i=0; i<len/16; i++) {

                ////////////////////////////////////////////////////////////////////////////////

                // calculate new sin, cos for next time through loop, using our incremental

                // updating algorithm

                ////////////////////////////////////////////////////////////////////////////////

                tempCos1    = newCos1 - (updateA*newCos1  + updateB * newSin1);

                newSin1     = newSin1 - (updateA*newSin1  - updateB * newCos1);

                newCos1     = tempCos1;

                tempCos2    = newCos2 - (updateA*newCos2  + updateB * newSin2);

                newSin2     = newSin2 - (updateA*newSin2  - updateB * newCos2);

                newCos2     = tempCos2;

                ////////////////////////////////////////////////////////////////////////////////

                // load input vectors

                ////////////////////////////////////////////////////////////////////////////////

                vLoBottom   = *pInLoBottom++;

                vHiBottom   = *pInHiBottom++;

                vLoTop      = *pInLoTop--;

                vHiTop      = *pInHiTop--;

                ////////////////////////////////////////////////////////////////////////////////

                // store new cos & sin to transition vector

                ////////////////////////////////////////////////////////////////////////////////

                ((float*)&vFTransition)[0] = newSin2;

                ((float*)&vFTransition)[1] = newSin1;

                ((float*)&vFTransition)[2] = newCos2;

                ((float*)&vFTransition)[3] = newCos1;

                vFNegTransition = vec_sub(vZero, vFTransition);

                ////////////////////////////////////////////////////////////////////////////////

                // calc sum of input data

                ////////////////////////////////////////////////////////////////////////////////

                *pOutLoBottom++     = vec_add(vLoBottom, vHiBottom);        

                *pOutLoTop--        = vec_add(vLoTop, vHiTop);      

                ////////////////////////////////////////////////////////////////////////////////

                // calc difference of input data

                ////////////////////////////////////////////////////////////////////////////////

                vDiffBottom     = vec_sub(vLoBottom, vHiBottom);

                vDiffTop        = vec_sub(vLoTop, vHiTop);

                ////////////////////////////////////////////////////////////////////////////////

                // multiply difference by twist factor

                ////////////////////////////////////////////////////////////////////////////////

                vTwistBottom    = vec_madd(vCosLo, vDiffBottom, vZero);     

                vTwistTop       = vec_madd(vCosHi, vDiffTop, vZero);        

                vSwappedDiffTop     = vec_perm(vDiffTop, vDiffTop, vSwapPairPerm);  

                vSwappedDiffBottom  = vec_perm(vDiffBottom, vDiffBottom, vSwapPairPerm);    

                vTwistBottom    = vec_madd(vSinLo, vSwappedDiffBottom, vTwistBottom);

                vTwistTop       = vec_madd(vSinHi, vSwappedDiffTop, vTwistTop);

                ////////////////////////////////////////////////////////////////////////////////

                // generate new cos & sin vectors from transition vector and previously

                // calculated steps.

                ////////////////////////////////////////////////////////////////////////////////

                vCosLo = vec_sub(vZero, vec_perm(vCosHi, vFNegTransition, vCosLoPerm));

                vCosHi = vec_mergel(vFNegTransition, vFNegTransition);

                vSinLo = vec_perm(vFTransition, vFNegTransition, vSinPerm);

                vSinLo = vec_sel(vSinHi, vSinLo, vHiPairSelect);

                vSinHi = vec_perm(vFTransition, vFNegTransition, vSinPerm);

                ////////////////////////////////////////////////////////////////////////////////

                // load input vectors for next sum & difference calculation

                ////////////////////////////////////////////////////////////////////////////////

                vLoBottom   = *pInLoBottom++;

                vHiBottom   = *pInHiBottom++;

                vLoTop      = *pInLoTop--;

                vHiTop      = *pInHiTop--;

                ////////////////////////////////////////////////////////////////////////////////

                // store twist-multiplied difference.

                ////////////////////////////////////////////////////////////////////////////////

                *pOutHiBottom++     = vTwistBottom;

                *pOutHiTop--        = vTwistTop;

                ////////////////////////////////////////////////////////////////////////////////

                // calculate new sin, cos for next time through loop, using our incremental

                // updating algorithm

                ////////////////////////////////////////////////////////////////////////////////

                tempCos1    = newCos1 - (updateA*newCos1  + updateB * newSin1);

                newSin1     = newSin1 - (updateA*newSin1  - updateB * newCos1);

                newCos1     = tempCos1;

                tempCos2    = newCos2 - (updateA*newCos2  + updateB * newSin2);

                newSin2     = newSin2 - (updateA*newSin2  - updateB * newCos2);

                newCos2     = tempCos2;

                ////////////////////////////////////////////////////////////////////////////////

                // store new cos & sin to transition vector

                ////////////////////////////////////////////////////////////////////////////////

                ((float*)&vFTransition)[0] = newSin2;

                ((float*)&vFTransition)[1] = newSin1;

                ((float*)&vFTransition)[2] = newCos2;

                ((float*)&vFTransition)[3] = newCos1;

                vFNegTransition = vec_sub(vZero, vFTransition);

                ////////////////////////////////////////////////////////////////////////////////

                // calc sum of input data

                ////////////////////////////////////////////////////////////////////////////////

                *pOutLoBottom++     = vec_add(vLoBottom, vHiBottom);        

                *pOutLoTop--        = vec_add(vLoTop, vHiTop);      

                ////////////////////////////////////////////////////////////////////////////////

                // calc difference of input data

                ////////////////////////////////////////////////////////////////////////////////

                vDiffBottom     = vec_sub(vLoBottom, vHiBottom);

                vDiffTop        = vec_sub(vLoTop, vHiTop);

                ////////////////////////////////////////////////////////////////////////////////

                // multiply difference by twist factor

                ////////////////////////////////////////////////////////////////////////////////

                vTwistBottom    = vec_madd(vCosLo, vDiffBottom, vZero);     

                vTwistTop       = vec_madd(vCosHi, vDiffTop, vZero);        

                vSwappedDiffTop     = vec_perm(vDiffTop, vDiffTop, vSwapPairPerm);  

                vSwappedDiffBottom  = vec_perm(vDiffBottom, vDiffBottom, vSwapPairPerm);    

                vTwistBottom    = vec_madd(vSinLo, vSwappedDiffBottom, vTwistBottom);

                vTwistTop       = vec_madd(vSinHi, vSwappedDiffTop, vTwistTop);

                ////////////////////////////////////////////////////////////////////////////////

                // generate new cos & sin vectors from transition vector and previously

                // calculated steps.

                ////////////////////////////////////////////////////////////////////////////////

                vCosLo = vec_sub(vZero, vec_perm(vCosHi, vFNegTransition, vCosLoPerm));

                vCosHi = vec_mergel(vFNegTransition, vFNegTransition);

                vSinLo = vec_perm(vFTransition, vFNegTransition, vSinPerm);

                vSinLo = vec_sel(vSinHi, vSinLo, vHiPairSelect);

                vSinHi = vec_perm(vFTransition, vFNegTransition, vSinPerm);

                ////////////////////////////////////////////////////////////////////////////////

                // store twist-multiplied difference.

                ////////////////////////////////////////////////////////////////////////////////

                *pOutHiBottom++     = vTwistBottom;

                *pOutHiTop--        = vTwistTop;

            }   

        } else {

            ////////////////////////////////////////////////////////////////////////////////

            // We are doing an inverse FFT, so we must adjust output by dividing by 2.

            ////////////////////////////////////////////////////////////////////////////////

            vOneHalf        = (vector float)(.5);

            for (i=0; i<len/16; i++) {

                ////////////////////////////////////////////////////////////////////////////////

                // calculate new sin, cos for next time through loop, using our incremental

                // updating algorithm

                ////////////////////////////////////////////////////////////////////////////////

                tempCos1    = newCos1 - (updateA*newCos1  + updateB * newSin1);

                newSin1     = newSin1 - (updateA*newSin1  - updateB * newCos1);

                newCos1     = tempCos1;

                tempCos2    = newCos2 - (updateA*newCos2  + updateB * newSin2);

                newSin2     = newSin2 - (updateA*newSin2  - updateB * newCos2);

                newCos2     = tempCos2;

                ////////////////////////////////////////////////////////////////////////////////

                // load input vectors

                ////////////////////////////////////////////////////////////////////////////////

                vLoBottom   = *pInLoBottom++;

                vHiBottom   = *pInHiBottom++;

                vLoTop      = *pInLoTop--;

                vHiTop      = *pInHiTop--;

                ////////////////////////////////////////////////////////////////////////////////

                // do inverse adjusting by dividing all elements by two. 

                ////////////////////////////////////////////////////////////////////////////////

                vLoBottom = vec_madd(vLoBottom, vOneHalf, vZero);

                vHiBottom = vec_madd(vHiBottom, vOneHalf, vZero);

                vLoTop = vec_madd(vLoTop, vOneHalf, vZero);

                vHiTop = vec_madd(vHiTop, vOneHalf, vZero);

                ////////////////////////////////////////////////////////////////////////////////

                // store new cos & sin to transition vector

                ////////////////////////////////////////////////////////////////////////////////

                ((float*)&vFTransition)[0] = newSin2;

                ((float*)&vFTransition)[1] = newSin1;

                ((float*)&vFTransition)[2] = newCos2;

                ((float*)&vFTransition)[3] = newCos1;

                vFNegTransition = vec_sub(vZero, vFTransition);

                ////////////////////////////////////////////////////////////////////////////////

                // calc sum of input data

                ////////////////////////////////////////////////////////////////////////////////

                *pOutLoBottom++     = vec_add(vLoBottom, vHiBottom);        

                *pOutLoTop--        = vec_add(vLoTop, vHiTop);      

                ////////////////////////////////////////////////////////////////////////////////

                // calc difference of input data

                ////////////////////////////////////////////////////////////////////////////////

                vDiffBottom     = vec_sub(vLoBottom, vHiBottom);

                vDiffTop        = vec_sub(vLoTop, vHiTop);

                ////////////////////////////////////////////////////////////////////////////////

                // multiply difference by twist factor

                ////////////////////////////////////////////////////////////////////////////////

                vTwistBottom    = vec_madd(vCosLo, vDiffBottom, vZero);     

                vTwistTop       = vec_madd(vCosHi, vDiffTop, vZero);        

                vSwappedDiffTop     = vec_perm(vDiffTop, vDiffTop, vSwapPairPerm);  

                vSwappedDiffBottom  = vec_perm(vDiffBottom, vDiffBottom, vSwapPairPerm);    

                vTwistBottom    = vec_madd(vSinLo, vSwappedDiffBottom, vTwistBottom);

                vTwistTop       = vec_madd(vSinHi, vSwappedDiffTop, vTwistTop);

                ////////////////////////////////////////////////////////////////////////////////

                // generate new cos & sin vectors from transition vector and previously

                // calculated steps.

                ////////////////////////////////////////////////////////////////////////////////

                vCosLo = vec_sub(vZero, vec_perm(vCosHi, vFNegTransition, vCosLoPerm));

                vCosHi = vec_mergel(vFNegTransition, vFNegTransition);

                vSinLo = vec_perm(vFTransition, vFNegTransition, vSinPerm);

                vSinLo = vec_sel(vSinHi, vSinLo, vHiPairSelect);

                vSinHi = vec_perm(vFTransition, vFNegTransition, vSinPerm);

                ////////////////////////////////////////////////////////////////////////////////

                // store twist-multiplied difference.

                ////////////////////////////////////////////////////////////////////////////////

                *pOutHiBottom++     = vTwistBottom;

                *pOutHiTop--        = vTwistTop;

                ////////////////////////////////////////////////////////////////////////////////

                // calculate new sin, cos for next time through loop, using our incremental

                // updating algorithm

                ////////////////////////////////////////////////////////////////////////////////

                tempCos1    = newCos1 - (updateA*newCos1  + updateB * newSin1);

                newSin1     = newSin1 - (updateA*newSin1  - updateB * newCos1);

                newCos1     = tempCos1;

                tempCos2    = newCos2 - (updateA*newCos2  + updateB * newSin2);

                newSin2     = newSin2 - (updateA*newSin2  - updateB * newCos2);

                newCos2     = tempCos2;

                ////////////////////////////////////////////////////////////////////////////////

                // load input vectors

                ////////////////////////////////////////////////////////////////////////////////

                vLoBottom   = *pInLoBottom++;

                vHiBottom   = *pInHiBottom++;

                vLoTop      = *pInLoTop--;

                vHiTop      = *pInHiTop--;

                ////////////////////////////////////////////////////////////////////////////////

                // do inverse adjusting by dividing all elements by two. 

                ////////////////////////////////////////////////////////////////////////////////

                vLoBottom = vec_madd(vLoBottom, vOneHalf, vZero);

                vHiBottom = vec_madd(vHiBottom, vOneHalf, vZero);

                vLoTop = vec_madd(vLoTop, vOneHalf, vZero);

                vHiTop = vec_madd(vHiTop, vOneHalf, vZero);

                ////////////////////////////////////////////////////////////////////////////////

                // store new cos & sin to transition vector

                ////////////////////////////////////////////////////////////////////////////////

                ((float*)&vFTransition)[0] = newSin2;

                ((float*)&vFTransition)[1] = newSin1;

                ((float*)&vFTransition)[2] = newCos2;

                ((float*)&vFTransition)[3] = newCos1;

                vFNegTransition = vec_sub(vZero, vFTransition);

                ////////////////////////////////////////////////////////////////////////////////

                // calc sum of input data

                ////////////////////////////////////////////////////////////////////////////////

                *pOutLoBottom++     = vec_add(vLoBottom, vHiBottom);        

                *pOutLoTop--        = vec_add(vLoTop, vHiTop);      

                ////////////////////////////////////////////////////////////////////////////////

                // calc difference of input data

                ////////////////////////////////////////////////////////////////////////////////

                vDiffBottom     = vec_sub(vLoBottom, vHiBottom);

                vDiffTop        = vec_sub(vLoTop, vHiTop);

                ////////////////////////////////////////////////////////////////////////////////

                // multiply difference by twist factor

                ////////////////////////////////////////////////////////////////////////////////

                vTwistBottom    = vec_madd(vCosLo, vDiffBottom, vZero);     

                vTwistTop       = vec_madd(vCosHi, vDiffTop, vZero);        

                vSwappedDiffTop     = vec_perm(vDiffTop, vDiffTop, vSwapPairPerm);  

                vSwappedDiffBottom  = vec_perm(vDiffBottom, vDiffBottom, vSwapPairPerm);    

                vTwistBottom    = vec_madd(vSinLo, vSwappedDiffBottom, vTwistBottom);

                vTwistTop       = vec_madd(vSinHi, vSwappedDiffTop, vTwistTop);

                ////////////////////////////////////////////////////////////////////////////////

                // generate new cos & sin vectors from transition vector and previously

                // calculated steps.

                ////////////////////////////////////////////////////////////////////////////////

                vCosLo = vec_sub(vZero, vec_perm(vCosHi, vFNegTransition, vCosLoPerm));

                vCosHi = vec_mergel(vFNegTransition, vFNegTransition);

                vSinLo = vec_perm(vFTransition, vFNegTransition, vSinPerm);

                vSinLo = vec_sel(vSinHi, vSinLo, vHiPairSelect);

                vSinHi = vec_perm(vFTransition, vFNegTransition, vSinPerm);

                ////////////////////////////////////////////////////////////////////////////////

                // store twist-multiplied difference.

                ////////////////////////////////////////////////////////////////////////////////

                *pOutHiBottom++     = vTwistBottom;

                *pOutHiTop--        = vTwistTop;

            }   

        }

        ////////////////////////////////////////////////////////////////////////////////

        // perform N/2-length FFT on high half of data

        ////////////////////////////////////////////////////////////////////////////////

        result = FFTComplex(pData+len, len/2, isign);

        if (result == noErr) {

            ////////////////////////////////////////////////////////////////////////////////

            // perform N/2-length FFT on low half of data

            ////////////////////////////////////////////////////////////////////////////////

            result = FFTComplex(pData, len/2, isign);

            if (result == noErr) {

                ////////////////////////////////////////////////////////////////////////////////

                // Move low half of data to temporary buffer, which will leave room in

                // main data space to merge elements of low and high data in to the original

                // buffer space.

                ////////////////////////////////////////////////////////////////////////////////

                {                               

                    vector float *pSource = (vector float*)pData;

                    vector float *pDest = *(vector float**)gTempBufferHandle;

                    vector float v1, v2, v3, v4;

                    for (i = 0; i < len/16; i++) {

                        v1 = *pSource++;

                        v2 = *pSource++;

                        v3 = *pSource++;

                        v4 = *pSource++;

                        *pDest++ = v1;

                        *pDest++ = v2;

                        *pDest++ = v3;

                        *pDest++ = v4;

                    }

                }

                ////////////////////////////////////////////////////////////////////////////////

                // Merge low and high halves of data back into original data buffer.  Elements

                // are merged so that, given high and low signals X1 and X2, resulting merged

                // signal is X1(0), X2(0), X1(1), X2(1), ... , X1(N/2-1), X2(N/2-1)

                ////////////////////////////////////////////////////////////////////////////////

                pInHiBottom = (vector float*)*gTempBufferHandle;

                pInHiTop = (vector float*)(pData+len);

                pOutLoBottom = (vector float*)pData;    

                for (i=0; i<len/8; i++) {

                    vector float    vInEven1, vInOdd1;

                    vector float    vInEven2, vInOdd2;

                    vector float    vOut1, vOut2, vOut3, vOut4;

                    vInEven1 = *pInHiBottom++;

                    vInOdd1 = *pInHiTop++;

                    vInEven2 = *pInHiBottom++;

                    vInOdd2 = *pInHiTop++;

                    vOut1 = vec_perm(vInEven1, vInOdd1, vMergeHiPairPerm);

                    vOut2 = vec_perm(vInEven1, vInOdd1, vMergeLoPairPerm);

                    vOut3 = vec_perm(vInEven2, vInOdd2, vMergeHiPairPerm);

                    vOut4 = vec_perm(vInEven2, vInOdd2, vMergeLoPairPerm);

                    *pOutLoBottom++ = vOut1;

                    *pOutLoBottom++ = vOut2;

                    *pOutLoBottom++ = vOut3;

                    *pOutLoBottom++ = vOut4;

                }

            }

        }

    }

    return result;

}

////////////////////////////////////////////////////////////////////////////////

//  fft_matrix_forward_columnwise

//

//  Performs a recursive 2D matrix FFT on the data pointed to by

// pData. It leaves data in columnwise order.

//

//  There are three steps to the 2D matrix FFT.  First, a FFT is performed on

// each column of the matrix.  Second, each element in the matrix is multiplied

// by a twist factor.  In a length N signal, element X(j, k) is multiplied by

// e^(+/- 2 pi i j k / N ). Finally, an FFT is performed on each row of the

// matrix.  The resulting data is the FFT of the original signal, stored in

// columnwise order.

//  

//

////////////////////////////////////////////////////////////////////////////////

static OSErr    fft_matrix_forward_columnwise(float *pData, long length)

{

    long        i, j;

    long        pow;

    float       *pRow;

    long        rowCount;

    long        colCount;

    long        rowPower;

    long        colPower;

    Handle      rowBufferHandle = nil;

    OSErr       result = noErr;

    pow = log2max(length);

    ////////////////////////////////////////////////////////////////////////////////

    // length must be an exact power of 2

    ////////////////////////////////////////////////////////////////////////////////

    if ((1 << pow) != length) {

        return paramErr;

    }

    ////////////////////////////////////////////////////////////////////////////////

    // data must be 16-byte aligned because of vector load/store requirements

    ////////////////////////////////////////////////////////////////////////////////

    if (((long)pData) & 0x0F) {

        return paramErr;

    }

    ////////////////////////////////////////////////////////////////////////////////

    // calculate dimensions of matrix.  If matrix can not be square, then the

    // column count will be 2X the row count    

    ////////////////////////////////////////////////////////////////////////////////

    rowPower = pow / 2;

    colPower = rowPower;

    if (pow & 1) colPower++;

    rowCount = 1<<rowPower;

    colCount = 1<<colPower;

    ////////////////////////////////////////////////////////////////////////////////

    // allocate a buffer that will be used for column ferrying, and can store

    // two columns at a time.

    ////////////////////////////////////////////////////////////////////////////////

    rowBufferHandle = NewHandle(2*2*sizeof(float)*rowCount);

    result = MemError();

    if (result == noErr) {

        vector float            *pCurrentColumn;

        vector float            vIn1, vIn2;

        vector float            vOut1, vOut2;

        vector unsigned char    vMergeHiPairPerm;

        vector unsigned char    vMergeLoPairPerm;

        vector float            *pSplit1, *pSplit2;

        float                   *pFFT1, *pFFT2;

        vector unsigned char    vSwappedPerm;

        vector float            vSinSignMultiplier;

        vector float            vCosTwistA0;

        vector float            vCosTwistA1;

        vector float            vCosTemp0;

        vector float            vCosTemp1;

        vector float            vSinTwistA0;

        vector float            vSinTwistA1;

        vector float            vA;

        vector float            vB;

        vector float            vTransition;

        vector float            vZero = (vector float)(0);

        double                  fTemp;

        double                  baseAngle1;

        double                  baseAngle2;

        double                  fSinBaseAngle1;

        double                  fSinBaseAngle2;

        ////////////////////////////////////////////////////////////////////////////////

        // initialize a permute vector that, given input vectors:

        //

        //  X = x0 x1 x2 x3

        //  Y = y0 y1 y2 y3

        //

        // will generate

        //  Z = x0 x1 y0 y1

        ////////////////////////////////////////////////////////////////////////////////

        vMergeHiPairPerm = (vector unsigned char)(0, 1, 2, 3, 4, 5, 6, 7, 16, 17, 18, 19, 20, 21, 22, 23);

        ////////////////////////////////////////////////////////////////////////////////

        // initialize a permute vector that, given input vectors:

        //

        //  X = x0 x1 x2 x3

        //  Y = y0 y1 y2 y3

        //

        // will generate

        //  Z = x2 x3 y2 y3

        ////////////////////////////////////////////////////////////////////////////////

        vMergeLoPairPerm = (vector unsigned char)(8, 9, 10, 11, 12, 13, 14, 15, 24, 25, 26, 27, 28, 29, 30, 31);

        ////////////////////////////////////////////////////////////////////////////////

        // initialize a permute vector that, given input vectors:

        //

        //  X = x0 x1 x2 x3

        //  Y = y0 y1 y2 y3

        //

        // will generate

        //  Z = x1 x0 x3 y2

        ////////////////////////////////////////////////////////////////////////////////

        vSwappedPerm  = (vector unsigned char)(4, 5, 6, 7, 0, 1, 2, 3, 12, 13, 14, 15, 8, 9, 10, 11);

        ////////////////////////////////////////////////////////////////////////////////

        // vSinSignMultiplier is a vector used to create a vector in the form

        // sin(x) -sin(x) sin(y) -sin(y)

        ////////////////////////////////////////////////////////////////////////////////

        vSinSignMultiplier = (vector float)(1, -1, 1, -1);  

        ////////////////////////////////////////////////////////////////////////////////

        // we now loop through all columns, two columns at a time, and perform column

        // ferrying FFTs on the columns.  We first copy the columns to a separate

        // buffer, then perform the FFTs, and finally copy them back to their original

        // column positions

        ////////////////////////////////////////////////////////////////////////////////

        for (i=0; i<colCount/2; i++) {

            ////////////////////////////////////////////////////////////////////////////////

            // set up vectors for incremental updating of cos & sin vectors

            //

            // Given c = cos(w) and s = sin(w), if we want to find the cos and sin of w plus

            // some small angle d, then we can do so by defining:

            //  

            // a = 2 * ((sin(d/2)) ^ 2)

            // b = sin(d)

            // 

            // Then, we can calculate the cos and sin of the updated angle w+d as:

            // 

            // cos(w+d) = c - ac - bs

            // sin(w+d) = s - as + bc

            // 

            // the vectors vA and vB contain these increment factors a and b.  Note that

            // each outer loop iteration requires different increment angles.  Looking at

            // the matrix, the angle multipliers that we use for our cos & sin twist

            // multipliers look something like:

            //

            //      0/N     0/N     0/N     0/N     0/N     0/N     0/N ... 0/N

            //      0/N     1/N     2/N     3/N     4/N     5/N     6/N     

            //      0/N     2/N     4/N     6/N     8/N     10/N    12/N    .

            //      0/N     3/N     6/N     9/N     12/N    15/N    18/N    .   

            //      0/N     4/N     8/N     12/N    16/N    20/N    24/N

            //      .       .

            //      .       .

            //      0/N                                         ...         (N-1)(N-1)

            //

            // Since we're doing two columns at a time, the vA and vB vectors

            // need to have two different sets of increment values.  For example, for

            // the third and fourth columns (handled in the second time through the

            // outer loop) each row increments by 2/N and 3/N.  Since we actually 

            // have two sets of twist multiplier vectors, we skip a row for each, so

            // the increment angle is doubled.  This helps reduce calculation dependencies

            // and cuts the number of increment iterations in half, which helps reduce

            // errors in the single-precision calculations used.

            //      

            ////////////////////////////////////////////////////////////////////////////////

            baseAngle1 = 4*i*PI/(colCount*rowCount);

            baseAngle2 = (2*(2*i+1)*PI)/(colCount*rowCount);

            fSinBaseAngle1 = sin(baseAngle1);

            fTemp = 2*fSinBaseAngle1*fSinBaseAngle1;

            ((float*)&vTransition)[0] = fTemp;

            ((float*)&vTransition)[2] = sin(2*baseAngle1);

            fSinBaseAngle2 = sin(baseAngle2);

            fTemp = 2*fSinBaseAngle2*fSinBaseAngle2;

            ((float*)&vTransition)[1] = fTemp;

            ((float*)&vTransition)[3] = sin(2*baseAngle2);

            vA = vec_mergeh(vTransition, vTransition);

            vB = vec_mergel(vTransition, vTransition);

            vB = vec_madd(vB, vSinSignMultiplier, vZero);

            ////////////////////////////////////////////////////////////////////////////////

            //  Set up the initial twist multiplier vectors for the beginning of this

            // column.  Values are calculated in the scalar domain, and then

            // transferred to the vector domain via the vTransition vector.

            ////////////////////////////////////////////////////////////////////////////////

            ((float*)&vTransition)[0] = 1;

            ((float*)&vTransition)[2] = cos(baseAngle1);

            ((float*)&vTransition)[3] = cos(baseAngle2);

            vCosTwistA0 = vec_splat(vTransition, 0);

            vCosTwistA1 = vec_mergel(vTransition, vTransition);

            vSinTwistA0 = vZero;

            ((float*)&vTransition)[2] = fSinBaseAngle1;

            ((float*)&vTransition)[3] = fSinBaseAngle2;

            vSinTwistA1 = vec_mergel(vTransition, vTransition);

            vSinTwistA1 = vec_madd(vSinTwistA1, vSinSignMultiplier, vZero);

            ////////////////////////////////////////////////////////////////////////////////

            // initialize pointers in temporary buffer for storing copied columns.

            ////////////////////////////////////////////////////////////////////////////////

            pSplit1 = *(vector float**)rowBufferHandle;

            pSplit2 = pSplit1 + rowCount / 2;

            ////////////////////////////////////////////////////////////////////////////////

            // point at top of columns to be copied

            ////////////////////////////////////////////////////////////////////////////////

            pCurrentColumn = ((vector float*)pData) + i;

            ////////////////////////////////////////////////////////////////////////////////

            // loop through all rows in the current column, copying elements to the

            // temporary buffer.  We load two rows at a time, and permute the vectors

            // to create a single vector that contains the appropriate column elements

            // from both rows.

            ////////////////////////////////////////////////////////////////////////////////

            for (j=0; j < rowCount/2; j++) {

                ////////////////////////////////////////////////////////////////////////////////

                // read in two rows worth of vectors (two rows X two columns)

                ////////////////////////////////////////////////////////////////////////////////

                vIn1 = *pCurrentColumn;

                pCurrentColumn += colCount/2;

                vIn2 = *pCurrentColumn;

                pCurrentColumn += colCount/2;

                ////////////////////////////////////////////////////////////////////////////////

                // permute the vectors so that the two left column entries are in one vector

                // and the two right column entries are in one vector

                ////////////////////////////////////////////////////////////////////////////////

                vOut1 = vec_perm(vIn1, vIn2, vMergeHiPairPerm);

                vOut2 = vec_perm(vIn1, vIn2, vMergeLoPairPerm);

                ////////////////////////////////////////////////////////////////////////////////

                // store the split columns to the appropriate buffers

                ////////////////////////////////////////////////////////////////////////////////

                *pSplit1++ = vOut1;

                *pSplit2++ = vOut2;

            }

            ////////////////////////////////////////////////////////////////////////////////

            // perform FFTs on the two columns of data that we have copied to the temp

            // buffer

            ////////////////////////////////////////////////////////////////////////////////

            pFFT1 = *(float **)rowBufferHandle;

            pFFT2 = pFFT1 + 2*rowCount;

            result = FFTComplex(pFFT1, rowCount, -1);

            if (result != noErr) break;

            result = FFTComplex(pFFT2, rowCount, -1);

            if (result != noErr) break;

            ////////////////////////////////////////////////////////////////////////////////

            // point at the beginning of the two copied column buffers, and at the top

            // of the columns in the matrix where we will store them back.

            ////////////////////////////////////////////////////////////////////////////////

            pSplit1 = *(vector float**)rowBufferHandle;

            pSplit2 = pSplit1 + rowCount / 2;

            pCurrentColumn = ((vector float*)pData) + i;

            ////////////////////////////////////////////////////////////////////////////////

            // loop through all of the column entries that are stored in the temp buffer,

            // and merge them to be stored back into the columns of the matrix.  At the

            // same time, we perform the twist multiply on the entries.

            ////////////////////////////////////////////////////////////////////////////////

            for (j=0; j < rowCount/2; j++) {

                vector float    vNew1, vNew2;

                ////////////////////////////////////////////////////////////////////////////////

                // get two vectors of column data

                ////////////////////////////////////////////////////////////////////////////////

                vIn1 = *pSplit1++;

                vIn2 = *pSplit2++;

                ////////////////////////////////////////////////////////////////////////////////

                // turn them into row vectors

                ////////////////////////////////////////////////////////////////////////////////

                vNew1 = vec_perm(vIn1, vIn2, vMergeHiPairPerm);

                vNew2 = vec_perm(vIn1, vIn2, vMergeLoPairPerm);

                ////////////////////////////////////////////////////////////////////////////////

                // do twist multiplication on both row vectors

                ////////////////////////////////////////////////////////////////////////////////

                vOut1 = vec_madd(vCosTwistA0, vNew1, vZero);

                vOut1 = vec_madd(vSinTwistA0, vec_perm(vNew1, vNew1, vSwappedPerm), vOut1);

                vOut2 = vec_madd(vCosTwistA1, vNew2, vZero);

                vOut2 = vec_madd(vSinTwistA1, vec_perm(vNew2, vNew2, vSwappedPerm), vOut2);

                ////////////////////////////////////////////////////////////////////////////////

                // store row vectors back into the matrix

                ////////////////////////////////////////////////////////////////////////////////

                *pCurrentColumn = vOut1;

                pCurrentColumn += colCount / 2;

                *pCurrentColumn = vOut2;                

                pCurrentColumn += colCount / 2;

                ////////////////////////////////////////////////////////////////////////////////

                // update sin & cos vectors for twist multiply

                ////////////////////////////////////////////////////////////////////////////////

                vCosTemp0 = vec_nmsub(vCosTwistA0, vA, vCosTwistA0);

                vCosTemp0 = vec_nmsub(vSinTwistA0, vB, vCosTemp0);

                vSinTwistA0 = vec_nmsub(vSinTwistA0, vA, vSinTwistA0);

                vSinTwistA0 = vec_madd (vCosTwistA0, vB, vSinTwistA0);

                vCosTwistA0 = vCosTemp0;

                vCosTemp1 = vec_nmsub(vCosTwistA1, vA, vCosTwistA1);

                vCosTemp1 = vec_nmsub(vSinTwistA1, vB, vCosTemp1);

                vSinTwistA1 = vec_nmsub(vSinTwistA1, vA, vSinTwistA1);

                vSinTwistA1 = vec_madd (vCosTwistA1, vB, vSinTwistA1);

                vCosTwistA1 = vCosTemp1;

            }

        }

        if (result == noErr) {

            ////////////////////////////////////////////////////////////////////////////////

            // finally, loop through all rows of matrix, performing an FFT on each row

            ////////////////////////////////////////////////////////////////////////////////

            pRow = pData;

            for (i=rowCount-1; i>=0; i--) {