c代码实现 ifft运算_FFT实现的C语言代码（转载）-（基2FFT及IFFT算法C语言实现）...

FFT实现的C语言代码- -(基2FFT及IFFT算法C语言实现)

Given two images A and B, use image B to cover image A. Where would we put B on A, so that the overlapping part of A and B has the most likelihood? To simplify the problem, we assume that A and B only contain numbers between 0 and 255. The difference between A and B is defined as the square sum of the differences of corresponding elements in the overlapped parts of A and B.

For example, we have

A (3 * 3):

a1

a2

a3

B (2 * 2):

b1

b2

a4

a5

a6

b4

b5

a7

a8

a9

When B is placed on position a5, the difference of them is ((b1-a5)^2 + (b2-a6)^2 + (b4-a8)^2 + (b5-a9)^2). Now we hope to have the position of the top left corner of B that gives the minimum difference. (B must completely reside on A)

It is clear that a simple solution will appear with very low efficiency when A and B have too many elements. But we can use 1-dimensional repeat convolution, which can be computed by Fast Fourier Transform (FFT), to improve the performance.

A program with explanation of FFT is given below:

/**

* Given two sequences {a1, a2, a3.. an} and {b1, b2, b3... bn},

* their repeat convolution means:

* r1 = a1*b1 + a2*b2 + a3*b3 + ... + an*bn

* r2 = a1*bn + a2*b1 + a3*b2 + ... + an*bn-1

* r3 = a1*bn-1 + a2*bn + a3*b1 + ... + an*bn-2

* ...

* rn = a1*b2 + a2*b3 + a3*b4 + ... + an-1*bn + an*b1

* Notice n >= 2 and n must be power of 2.

*/

#include

#include

#include

#define for if (0); else for

using namespace std;

const int MaxFastBits = 16;

int **gFFTBitTable = 0;

int NumberOfBitsNeeded(int PowerOfTwo) {

for (int i = 0;; ++i) {

if (PowerOfTwo & (1 << i)) {

return i;

}

}

}

int ReverseBits(int index, int NumBits) {

int ret = 0;

for (int i = 0; i < NumBits; ++i, index >>= 1) {

ret = (ret << 1) | (index & 1);

}

return ret;

}

void InitFFT() {

gFFTBitTable = new int *[MaxFastBits];

for (int i = 1, length = 2; i <= MaxFastBits; ++i, length <<= 1) {

gFFTBitTable[i - 1] = new int[length];

for (int j = 0; j < length; ++j) {

gFFTBitTable[i - 1][j] = ReverseBits(j, i);

}

}

}

inline int FastReverseBits(int i, int NumBits) {

return NumBits <= MaxFastBits ? gFFTBitTable[NumBits - 1][i] : ReverseBits(i, NumBits);

}

void FFT(bool InverseTransform, vector >& In, vector >& Out) {

if (!gFFTBitTable) { InitFFT(); }

// simultaneous data copy and bit-reversal ordering into outputs

int NumSamples = In.size();

int NumBits = NumberOfBitsNeeded(NumSamples);

for (int i = 0; i < NumSamples; ++i) {

Out[FastReverseBits(i, NumBits)] = In[i];

}

// the FFT process

double angle_numerator = acos(-1.) * (InverseTransform ? -2 : 2);

for (int BlockEnd = 1, BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {

double delta_angle = angle_numerator / BlockSize;

double sin1 = sin(-delta_angle);

double cos1 = cos(-delta_angle);

double sin2 = sin(-delta_angle * 2);

double cos2 = cos(-delta_angle * 2);

for (int i = 0; i < NumSamples; i += BlockSize) {

complex a1(cos1, sin1), a2(cos2, sin2);

for (int j = i, n = 0; n < BlockEnd; ++j, ++n) {

complex a0(2 * cos1 * a1.real() - a2.real(), 2 * cos1 * a1.imag() - a2.imag());

a2 = a1;

a1 = a0;

complex a = a0 * Out[j + BlockEnd];

Out[j + BlockEnd] = Out[j] - a;

Out[j] += a;

}

}

BlockEnd = BlockSize;

}

// normalize if inverse transform

if (InverseTransform) {

for (int i = 0; i < NumSamples; ++i) {

Out[i] /= NumSamples;

}

}

}

vector convolution(vector a, vector b) {

int n = a.size();

vector > s(n), d1(n), d2(n), y(n);

for (int i = 0; i < n; ++i) {

s[i] = complex(a[i], 0);

}

FFT(false, s, d1);

s[0] = complex(b[0], 0);

for (int i = 1; i < n; ++i) {

s[i] = complex(b[n - i], 0);

}

FFT(false, s, d2);

for (int i = 0; i < n; ++i) {

y[i] = d1[i] * d2[i];

}

FFT(true, y, s);

vector ret(n);

for (int i = 0; i < n; ++i) {

ret[i] = s[i].real();

}

return ret;

}

int main() {

double a[4] = {1, 2, 3, 4}, b[4] = {1, 2, 3, 4};

vector r = convolution(vector(a, a + 4), vector(b, b + 4));

// r[0] = 30 (1*1 + 2*2 + 3*3 + 4*4)

// r[1] = 24 (1*4 + 2*1 + 3*2 + 4*3)

// r[2] = 22 (1*3 + 2*4 + 3*1 + 4*2)

// r[3] = 24 (1*2 + 2*3 + 3*4 + 4*1)

return 0;

}

Input

The first line contains n (1 <= n <= 10), the number of test cases.

For each test case, the first line contains four integers m, n, p and q, where A is a matrix of m * n, B is a matrix of p * q (2 <= m, n, p, q <= 500, m >= p, n >= q). The following m lines are the elements of A and p lines are the elements of B.

Output

For each case, print the position that gives the minimum difference (the top left corner of A is (1, 1)). You can assume that each test case has a unique solution.

Sample Input

2

2 2 2 2

1 2

3 4

2 3

1 4

3 3 2 2

0 5 5

0 5 5

0 0 0

5 5

5 5

Sample Output

1 1

1 2

Author: DU, Peng