c++ DFT 二维傅里叶变换

最近写了一个计算图像二维傅里叶变换的程序， c++， 采用索引表的方法计算。

程序如下

#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
using namespace cv;

#include <stdlib.h>
#include <stdio.h>

#define PI 3.14159265359

float *trigalesin;
float *trigalecos;

void calculate(int col, int row)
{
	trigalesin  = new float[col*row];
	trigalecos = new float[col*row];

	for (int y = 0; y< row; y++)
	{
		for(int x = 0; x< col; x++)
		{
			float tmp = -2*PI*(1.0*x/col+1.0*y/row);

			trigalesin[y*col+x] = sin(tmp);
			trigalecos[y*col+x] = cos(tmp);
		}
	}
}

void getDFT(int *a, int arows, int acols, float **b, float **c)
{
	

	*b = new float [arows*acols];
	*c = new float [arows*acols];

	for (int v = 0; v< arows; v++)
	{
		if (v*100%arows == 0)
		{
			printf("%d\n",v*100/arows );
		}
		for (int u = 0; u< acols; u++)
		{
			float tmpb = 0;
			float tmpc = 0;



			for (int y = 0; y< arows; y++)
			{
				for (int x = 0; x< acols; x++)
				{
					float t1 = a[y*acols+x];

					if (t1 == 0)
					{
						continue;
					}

					// float t2 = -2*PI*(u*x/acols+v*y/arows);

					int index1 = (u*x)%acols;
					int index2 = (v*y)%arows;
					int index = index2*acols+index1;


					tmpb += t1*trigalesin[index];
					tmpc += t1*trigalecos[index];
				}
			}

			(*b)[v*acols+u] = tmpb;
			(*c)[v*acols+u] = tmpc;
		}
	}
}

void getIDFT(float *a, float *b, float **c, int arows, int acols)
{
	*c = new float [arows*acols];

	for (int y = 0; y< arows; y++)
	{
		if (y*100%arows == 0)
		{
			printf("%d\n",y*100/arows );
		}

		for (int x = 0; x< acols; x++)
		{
			float tmpa = 0;
			float tmpb = 0;

			for (int v = 0; v< arows; v++)
			{
				for (int u = 0; u< acols; u++)
				{
					float tmpfa = a[v*acols+u];
					float tmpfb = b[v*acols+u];

					int index1 = (u*x)%acols;
					int index2 = (v*y)%arows;
					int index = index2*acols+index1;

					float tmpea = -trigalesin[index];//sin(2*PI*(u*x/acols+v*y/arows));
					float tmpeb = trigalecos[index];//cos(2*PI*(u*x/acols+v*y/arows));

					tmpa += tmpfa*tmpea - tmpfb*tmpeb;
					tmpb += tmpfa*tmpeb+tmpfb*tmpea;
				}
			}

			(*c)[y*acols+x] = sqrt(tmpa*tmpa+tmpb*tmpb);
		}
	}
}

void getPadding(Mat a, Mat b, int **adata, int **bdata, int &rows, int &cols)
{
	rows = pow(2, 1.0*(int)(log(b.rows*1.0)/log(2.0)+1));
	cols = pow(2, 1.0*(int)(log(b.cols*1.0)/log(2.0)+1));

	rows = b.rows;
	cols = b.cols;

	*adata = new int[rows*cols];
	*bdata =  new int[rows*cols];



	for (int i = 0; i< rows; i++)
	{
		for (int j = 0; j< cols; j++)
		{
			(*adata)[i*cols+j] = 0;
			(*bdata)[i*cols+j] = 0;
		}
	}

	for (int i = 0; i< a.rows; i++)
	{
		for (int j = 0; j< a.cols; j++)
		{
			(*adata)[i*cols+j] = a.at<uchar>(i, j);
		}
	}

	for (int i = 0; i< b.rows; i++)
	{
		for (int j = 0; j< b.cols; j++)
		{
			(*bdata)[i*cols+j] = b.at<uchar>(i, j);
		}
	}
}

void getMul(float *ab, float *ac, float *bb, float *bc, float **cb, float **cc, int rows, int cols)
{
	*cb = new float [rows*cols];
	*cc = new float [rows*cols];

	for (int i = 0; i< rows; i++)
	{
		for (int j = 0; j< cols; j++)
		{
			(*cb)[i*cols+j] = ab[i*cols+j]*bb[i*cols+j] - ac[i*cols+j]*bc[i*cols+j];
			(*cc)[i*cols+j] = ab[i*cols+j]*bc[i*cols+j]+ac[i*cols+j]*bb[i*cols+j];
		}
	}
}

over.

另外， dft做出来了， 但是如何计算图像位置还不清楚。