2

I am using iOS Accelerate framework for finding FFT of a 2D array. The code below works correctly only for power of 2 images. We have to pad input arrays with zeros for non power of 2 images. But I am not able to do the padding correctly. Currently I pad arrays like below 

float inputImg2D[3][3] = { 1,1,1, 1,1,1, 1,1,1 }; 
float paddedImg2D[4][4] = { 1,1,1,0, 1,1,1,0, 1,1,1,0, 0,0,0,0 };
float expectedOutput[6]6] = { 9,0,0,0,0,0
                              0,0,0,0,0,0
                              0,0,0,0,0,0
                              0,0,0,0,0,0
                              0,0,0,0,0,0
                              0,0,0,0,0,0 };

For a 4*4 array, I am correctly get output as 8*8 array with value 16 at (0,0).

Accelerate FFT Code. 

/* 
 * 2D fft sample working only for power of 2 images.
 * expected output for below 3*3 array is a 6*6 array with value 9 at ( 0,0) - all other values will be zero
 * expected output for below 4*4 array is a 8*8 array with value 16 at (0,0) - all other values will be zero
 */

#include <stdio.h>
#include "Accelerate/Accelerate.h"

#define NON_POWER_OF_2_TEST_WILL_FAIL

int main(int argc, const char * argv[]) {

#ifdef NON_POWER_OF_2_TEST_WILL_FAIL
    const int IMG_ROWS = 3;
    const int IMG_COLS = 3;
    float img2D[3][3] = { 1,1,1, 1,1,1, 1,1,1 };
#else
    const int IMG_ROWS = 4;
    const int IMG_COLS = 4;
    float img2D[4][4] = { 1,1,1,1, 1,1,1,1, 1,1,1,1, 1,1,1,1 };
#endif

    /*  build necessary values for fft setup */
    int maxDimension = ( IMG_ROWS > IMG_COLS ) ? IMG_ROWS : IMG_COLS;
    int optimalDftSize = ceil( log2( maxDimension) );

    /* fft setup */
    FFTSetup fftSetup = vDSP_create_fftsetup( optimalDftSize, FFT_RADIX2 );

    /* expand images to power of two size with zero values*/
    COMPLEX_SPLIT in_fft;
    int optimalDftWidth = 1 << optimalDftSize;
    int optimalDftHeight = 1 << optimalDftSize;
    int numElements = optimalDftWidth * optimalDftHeight;
    in_fft.realp = ( float* ) calloc ( numElements, sizeof(float) );
    in_fft.imagp = ( float* ) calloc ( numElements, sizeof(float) );

    /* assign image pixels if only in range */
    for ( int i = 0; i < optimalDftWidth; i++ ) {
        for ( int j = 0; j < optimalDftHeight; j++ ) {
            if (i < IMG_ROWS && j < IMG_COLS) {
                in_fft.realp[i * optimalDftHeight + j] = img2D[i][j];
                //in_fft.imagp[i] = 0.0;
            }
        }
    }

    /* do fft in place */
    int rowStride = 1;
    int columnStride = 0;
    vDSP_fft2d_zip(fftSetup, &in_fft, rowStride, columnStride, optimalDftSize, optimalDftSize, FFT_FORWARD);

    /* print results */
    for(int i=0; i < optimalDftWidth; ++i) {
        for(int j=0; j < optimalDftHeight; ++j) {
            printf (" %.2f, %.2f, ", in_fft.realp[i*optimalDftHeight+j], in_fft.imagp[i*optimalDftHeight+j] );
        }
        printf("\n");
    }

    /* TODO: free resources */
    return 0;
}
kiranpradeep
  • 10,859
  • 4
  • 50
  • 82
  • I think you're misunderstanding how zero padding works - the zeroes aren't ignored, they effectively produce interpolation in the frequency domain. So your code is (probably) working correctly - you're just misinterpreting the results. – Paul R Apr 02 '15 at 15:38
  • @PaulR Thanks. But, then how to get the `expectedOutput` ( as in question). For a 3*3 input array of all 1's, I am expecting a result of all zeros except a value of 9 at (0,0) in result array. Is my expectation wrong ? I assume the function `vDSP_fft2d_zip` expects padded zeros. – kiranpradeep Apr 02 '15 at 16:21
  • 1
    You'll only get the expected output if you use an actual 3x3 FFT. `vDSP_fft2d_zip` doesn't know that you're giving it padded data - it just blindly processes the data as a 4x4 FFT, so you get an interpolated 4x4 output. You'll need to use an FFT that supports non-power-of-2 sizes such as 3x3 (e.g. FFTW) or re-think whatever it is that you are trying to achieve with this operation. – Paul R Apr 02 '15 at 16:25
  • 1
    @PaulR Thanks. I assume, that meant, we cannot do an exact `n*n` dft with Accelerate framework, where `n` cannot be expressed as power of 2,3 or 5. The second part about rethinking also helped, since I was planning to do correlations with Accelerate which didn't really need the exact n*m dft. I will follow the `convolveDFT` algorithm as in http://docs.opencv.org/modules/core/doc/operations_on_arrays.html#dft. This comment has helped and let me know if could write it as an answer to be accepted. Or else I could do the same. – kiranpradeep Apr 03 '15 at 17:38
  • Glad to have helped - feel free to write up your conclusions as an answer for the benefit of future visitors and I'll upvote it. – Paul R Apr 03 '15 at 18:23

2 Answers2

3

You are likely doing the zero-padding correctly. Any expectation that the results will be the same after zero-padding is incorrect. All the FFT bin frequencies will instead be related to the larger sized array.

hotpaw2
  • 70,107
  • 14
  • 90
  • 153
  • 1
    Thanks. Yes. my expectation was wrong. Before posting, I was wondering, about how ios Accelerate was going to differentiate between the zero padded 4*4 matrix and the real 4*4 matrix of same data. Thought there could be some function parameters( possible strides ) to pass in `vDSP_fft2d_zip` to express the intent of 3*3 dft. It appears there is none. – kiranpradeep Apr 03 '15 at 17:36
1

Adding my assumptions made up from Paul R's comment responses to questions ( could be wrong as my understanding of FFT usage/purpose is nearly zero )

  1. With iOS Accelerate framework, you cannot exactly compute FFT of 2D array/matrix whose sizes are non powers of 2. You can zero pad the array to a power of 2, but that effects results.
  2. Mostly you could do with out non power of 2. For e.g I was trying to calculate convolution between two 2D matrices, and it wasn't necessary to do to get an exact m * n FFT. Sample here(1).
  3. Non related to the question: If you are trying to replace OpenCV's dft method with Accelerate vDSP_fft2d_zripD, in the hope of better performance, I feel there is none or only far worse. Even accepting reuse of buffer's are critical to Accelerate's performance, it was disappointing. May be my bad code.
kiranpradeep
  • 10,859
  • 4
  • 50
  • 82