This is more a "post a suggestion" topic rather than a question. And thank you if you are willing to read this whole.

I've been studing the code in the Nvidia Cuda SDK regarding how to operate a convolution between matrix and kernel using FFT (so basically you transform matrix and kernel and perform a multiplication in the frequency domain.. eventually you inverse-transform the result and you're done).

For the ones of you interested in the code quickly these are the classes http://nopaste.info/30c13e44fe.html http://nopaste.info/78d22afac2.html

Anyway I got everything that the code does... excepted one thing: the code seems to perform a pre-processing step before the multiplication. After the multiplication the code is again post-processed. So seems like a signal modulation or signal processing of some kind I can't understand.. if it was for me I'd perform the multiplication without all these complications, but I suppose the author of the code had good reasons to add these processing steps.

The routines are here:

inline __device__ void spPostprocessC2C(fComplex& D1, fComplex& D2, const fComplex& twiddle){
    float A1 = 0.5f * (D1.x + D2.x);
    float B1 = 0.5f * (D1.y - D2.y);
    float A2 = 0.5f * (D1.y + D2.y);
    float B2 = 0.5f * (D1.x - D2.x);

    D1.x = A1 + (A2 * twiddle.x + B2 * twiddle.y);
    D1.y = (A2 * twiddle.y - B2 * twiddle.x) + B1;
    D2.x = A1 - (A2 * twiddle.x + B2 * twiddle.y);
    D2.y = (A2 * twiddle.y - B2 * twiddle.x) - B1;

//Premultiply by 2 to account for 1.0 / (DZ * DY * DX) normalization
inline __device__ void spPreprocessC2C(fComplex& D1, fComplex& D2, const fComplex& twiddle){
    float A1 = /* 0.5f * */ (D1.x + D2.x);
    float B1 = /* 0.5f * */ (D1.y - D2.y);
    float A2 = /* 0.5f * */ (D1.y + D2.y);
    float B2 = /* 0.5f * */ (D1.x - D2.x);

    D1.x = A1 - (A2 * twiddle.x - B2 * twiddle.y);
    D1.y = (B2 * twiddle.x + A2 * twiddle.y) + B1;
    D2.x = A1 + (A2 * twiddle.x - B2 * twiddle.y);
    D2.y = (B2 * twiddle.x + A2 * twiddle.y) - B1;

inline __device__ void getTwiddle(fComplex& twiddle, float phase){
    __sincosf(phase, &twiddle.y, &twiddle.x);

and before calling these,

fComplex twiddle;
getTwiddle(twiddle, phaseBase * (float)x);
spPostprocessC2C(D1, D2, twiddle);

this code gets executed. x is the row index of the padded matrix, phaseBase is defined as follows

const double phaseBase = PI / (double)DX;

so it's PI / image_width. Seems like the code is "spreading" the twiddle factor for every row on the multiplication. I simply can't understand this need.. but I'm not an expert in signal processing.

Does anyone of you have any idea of the reason for this code to be run?

Thank you in advance for every comment/suggestion

  • $\begingroup$ If you have analyzed the code in detail, could you describe what it does in mathematical terms, or produce a non-CUDA version? Although it may be obvious to signal processing experts, I'm already lost at how this thing is called. I guess spProcess2D is the main function, and it operates on frequency-domain signals. Is this correct? Maybe the twiddling is just the last step of the input FFT and the first step of the output FFT, and was included in the CUDA kernel for performance reasons. $\endgroup$ – Sebastian Reichelt May 20 '11 at 21:37
  • $\begingroup$ Apparently, this document describes the implementation: developer.download.nvidia.com/compute/cuda/2_2/sdk/website/…. It doesn't seem to answer your question, though. $\endgroup$ – Sebastian Reichelt May 20 '11 at 21:43
  • $\begingroup$ Yes I've already seen that, without much success. The code does not implement the FFT because a library called "CUFFT" is called before the entire code I've posted. So the Cooley Tukey algorithm with prime general factorization is executed (and twiddles have got rid of too). The FFT code will be called once again to inverse the multiplication result, but the code I've posted comes in the scene before this step. In the code I've linked seems like there is a multiplication running, but the whole thing is just weird $\endgroup$ – Marco A. May 20 '11 at 22:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.