Why it is left to user to divide each element of an array by total number of points after computing backward FFT?

[audrius]

Hello,

Recently I 'discovered' that just after computing "DftiComputeForward(Desc_Handle, x_in)" and "DftiComputeBackward(Desc_Handle, x_in)" the array x_in is not equal to the initial: The user must divide each x_in element by total number of points in order to achieve equality.

Q1: why it is left to user to do that? note, that the user might not do that in the most efficient way thus reducing overall performance.

Q2: does user need to do something similar after computing only "DftiComputeForward(Desc_Handle, x_in)"?

Program code (VS2008 C++ Win32 MKL 10.0.5.025):
// 96.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include
#include "mkl_dfti.h"
#include
#include

int _tmain(int argc, _TCHAR* argv[])
{
printf("Hello World \n");

DFTI_DESCRIPTOR_HANDLE Desc_Handle = 0;
int Status, i1, i2, i3, m, i;
double d1, d2, d3, d4;
char *c;
double *x_in;

m = 128;
i1 = 1;

c = new char;
x_in = new double[2 * m];
i2 = 0;
i3 = m / 2;
for (i = 0; i < m; i++) { //Fills x_in array with gaussian pulse
d1 = 3.1622776601683795 * exp(-(i - i3) * (i - i3) * 0.013862943611198901);
x_in[i2] = d1;
x_in[i2 + 1] = 0.0;
i2 += 2;
}
i2 = 0;
d4 = 0;
for (i = 0; i < m; i++) { //searching for the peak value
d1 = x_in[i2];
d2 = x_in[i2 + 1];
d3 = d1 * d1 + d2 * d2;
if (d3 > d4) d4 = d3;
i2 += 2;
}
printf("%f \n", d4);
Status = DftiCreateDescriptor(&Desc_Handle, DFTI_DOUBLE, DFTI_COMPLEX, i1, m);
printf("DftiCreateDescriptor error_message = %s \n", DftiErrorMessage(Status));
Status = DftiCommitDescriptor( Desc_Handle );
printf("DftiCommitDescriptor error_message = %s \n", DftiErrorMessage(Status));
Status = DftiComputeForward(Desc_Handle, x_in);
printf("DftiComputeForward error_message = %s \n", DftiErrorMessage(Status));
Status = DftiComputeBackward(Desc_Handle, x_in);
printf("DftiComputeBackward error_message = %s \n", DftiErrorMessage(Status));
DftiFreeDescriptor(&Desc_Handle);
printf("DftiFreeDescriptor error_message = %s \n", DftiErrorMessage(Status));
i2 = 0;
for (i = 0; i < m; i++) { //dividing by number of points
x_in[i2] /= m;
x_in[i2+1] /= m;
i2 += 2;
}
i2 = 0;
d4 = 0;
for (i = 0; i < m; i++) { //searching for the peak value
d1 = x_in[i2];
d2 = x_in[i2 + 1];
d3 = d1 * d1 + d2 * d2;
if (d3 > d4) d4 = d3;
i2 += 2;
}
printf("%f \n", d4);
c = new char;
scanf(c);
return 0;
}

Thank you,

Audrius

[Vladimir_Petrov__Int]

Quoting - audrius

Hello,

Recently I 'discovered' that just after computing "DftiComputeForward(Desc_Handle, x_in)" and "DftiComputeBackward(Desc_Handle, x_in)" the array x_in is not equal to the initial: The user must divide each x_in element by total number of points in order to achieve equality.

Q1: why it is left to user to do that? note, that the user might not do that in the most efficient way thus reducing overall performance.

Q2: does user need to do something similar after computing only "DftiComputeForward(Desc_Handle, x_in)"?

Program code (VS2008 C++ Win32 MKL 10.0.5.025):
// 96.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include
#include "mkl_dfti.h"
#include
#include

int _tmain(int argc, _TCHAR* argv[])
{
printf("Hello World n");

DFTI_DESCRIPTOR_HANDLE Desc_Handle = 0;
int Status, i1, i2, i3, m, i;
double d1, d2, d3, d4;
char *c;
double *x_in;

m = 128;
i1 = 1;

c = new char;
x_in = new double[2 * m];
i2 = 0;
i3 = m / 2;
for (i = 0; i < m; i++) { //Fills x_in array with gaussian pulse
d1 = 3.1622776601683795 * exp(-(i - i3) * (i - i3) * 0.013862943611198901);
x_in[i2] = d1;
x_in[i2 + 1] = 0.0;
i2 += 2;
}
i2 = 0;
d4 = 0;
for (i = 0; i < m; i++) { //searching for the peak value
d1 = x_in[i2];
d2 = x_in[i2 + 1];
d3 = d1 * d1 + d2 * d2;
if (d3 > d4) d4 = d3;
i2 += 2;
}
printf("%f n", d4);
Status = DftiCreateDescriptor(&Desc_Handle, DFTI_DOUBLE, DFTI_COMPLEX, i1, m);
printf("DftiCreateDescriptor error_message = %s n", DftiErrorMessage(Status));
Status = DftiCommitDescriptor( Desc_Handle );
printf("DftiCommitDescriptor error_message = %s n", DftiErrorMessage(Status));
Status = DftiComputeForward(Desc_Handle, x_in);
printf("DftiComputeForward error_message = %s n", DftiErrorMessage(Status));
Status = DftiComputeBackward(Desc_Handle, x_in);
printf("DftiComputeBackward error_message = %s n", DftiErrorMessage(Status));
DftiFreeDescriptor(&Desc_Handle);
printf("DftiFreeDescriptor error_message = %s n", DftiErrorMessage(Status));
i2 = 0;
for (i = 0; i < m; i++) { //dividing by number of points
x_in[i2] /= m;
x_in[i2+1] /= m;
i2 += 2;
}
i2 = 0;
d4 = 0;
for (i = 0; i < m; i++) { //searching for the peak value
d1 = x_in[i2];
d2 = x_in[i2 + 1];
d3 = d1 * d1 + d2 * d2;
if (d3 > d4) d4 = d3;
i2 += 2;
}
printf("%f n", d4);
c = new char;
scanf(c);
return 0;
}

Thank you,

Audrius

Audrius,

In order to achieve the behaviour described in Q1 you need to set the DFTI_BACKWARD_SCALE parameter before committing the descriptor like this:

Status = DftiSetValue(Desc_Handle, DFTI_BACKWARD_SCALE, 1.0/(double)m);

Please note also that the DFTI interface uses MKL_LONG as the type of the return values as well as the integer parameters (most important when you pass integer arrays to DFTI).

Regarding Q2, could you be more specific as to what is "something similar"?

Best regards,

-Vladimir