怎样做二维付利叶变换？

Question

我在做雷达成像是要用到二维付利叶变换，未曾学过，不知该如何入手。请大家指点一下
输入信号为二维矩阵，是雷达在不同角度的回波信号，第一维是某角度的回波信号，第二维是转动角度值，如何对该矩阵做付利叶变换

Answer 1

void FreqChange(complex<double> * FD,long lWidth,long lHeight)
{//由于经过FFT变换时候，低频不在中心，变换 long i,j;
complex<double> temp;//移位临时变量

for(i = 0; i < lHeight/2; i++)
{
// 列
for(j = 0; j < lWidth; j++)
{
temp = FD[ (lHeight-1-i)*lWidth + j]; FD[ (lHeight-1-i)*lWidth + j] 
= FD[ (lHeight/2-1-i)*lWidth + (lWidth/2+j)%(lWidth)];

FD[ (lHeight/2-1-i)*lWidth + (lWidth/2+j)%(lWidth)]
= temp;
}
}}
/*************************************************************************
 *
 * 函数名称：
 *   FFT()
 *
 * 参数:
 *   complex<double> * TD - 指向时域数组的指针
 *   complex<double> * FD - 指向频域数组的指针
 *   r －2的幂数，即迭代次数
 *
 * 返回值:
 *   无。
 *
 * 说明:
 *   该函数用来实现快速付立叶变换。
 *
 ************************************************************************/void FFT(complex<double> * TD, complex<double> * FD, int r)
{
// 付立叶变换点数
long count;

// 循环变量
int i,j,k;

// 中间变量
int bfsize,p;

// 角度
double angle;

complex<double> *W,*X1,*X2,*X;

// 计算付立叶变换点数
count = 1 << r;

// 分配运算所需存储器
W  = new complex<double>[count / 2];
X1 = new complex<double>[count];
X2 = new complex<double>[count];

// 计算加权系数
for(i = 0; i < count / 2; i++)
{
angle = -i * PI * 2 / count;
W[i] = complex<double> (cos(angle), sin(angle));
}

// 将时域点写入X1
memcpy(X1, TD, sizeof(complex<double>) * count);

// 采用蝶形算法进行快速付立叶变换
for(k = 0; k < r; k++)
{
for(j = 0; j < 1 << k; j++)
{
bfsize = 1 << (r-k);
for(i = 0; i < bfsize / 2; i++)
{
p = j * bfsize;
X2[i + p] = X1[i + p] + X1[i + p + bfsize / 2];
X2[i + p + bfsize / 2] = (X1[i + p] - X1[i + p + bfsize / 2]) * W[i * (1<<k)];
}
}
X  = X1;
X1 = X2;
X2 = X;
}

// 重新排序
for(j = 0; j < count; j++)
{
p = 0;
for(i = 0; i < r; i++)
{
if (j&(1<<i))
{
p+=1<<(r-i-1);
}
}
FD[j]=X1[p];
}

// 释放内存
delete W;
delete X1;
delete X2;
}/*************************************************************************
 *
 * 函数名称：
 *   IFFT()
 *
 * 参数:
 *   complex<double> * FD - 指向频域值的指针
 *   complex<double> * TD - 指向时域值的指针
 *   r －2的幂数
 *
 * 返回值:
 *   无。
 *
 * 说明:
 *   该函数用来实现快速付立叶反变换。
 *
 ************************************************************************/void IFFT(complex<double> * FD, complex<double> * TD, int r)
{
// 付立叶变换点数
long count;

// 循环变量
int i;

complex<double> *X;

// 计算付立叶变换点数
count = 1 << r;

// 分配运算所需存储器
X = new complex<double>[count];

// 将频域点写入X
memcpy(X, FD, sizeof(complex<double>) * count);

// 求共轭
for(i = 0; i < count; i++)
{
X[i] = complex<double> (X[i].real(), -X[i].imag());
}

// 调用快速付立叶变换
FFT(X, TD, r);

// 求时域点的共轭
for(i = 0; i < count; i++)
{
TD[i] = complex<double> (TD[i].real() / count, -TD[i].imag() / count);
}

// 释放内存
delete X;
}/*************************************************************************
 *
 * 函数名称：
 *   DCT()
 *
 * 参数:
 *   double * f - 指向时域值的指针
 *   double * F - 指向频域值的指针
 *   r －2的幂数
 *
 * 返回值:
 *   无。
 *
 * 说明:
 *   该函数用来实现快速离散余弦变换。该函数利用2N点的快速付立叶变换
 * 来实现离散余弦变换。
 *
 ************************************************************************/
void  DCT(double *f, double *F, int r)
{
// 离散余弦变换点数
long count;

// 循环变量
int i;

// 中间变量
double dTemp;

complex<double> *X;

// 计算离散余弦变换点数
count = 1<<r;

// 分配内存
X = new complex<double>[count*2];

// 赋初值为0
memset(X, 0, sizeof(complex<double>) * count * 2);

// 将时域点写入数组X
for(i=0;i<count;i++)
{
X[i] = complex<double> (f[i], 0);
}

// 调用快速付立叶变换
FFT(X,X,r+1);

// 调整系数
dTemp = 1/sqrt(count);

// 求F[0]
F[0] = X[0].real() * dTemp;

dTemp *= sqrt(2);

// 求F[u]
for(i = 1; i < count; i++)
{
F[i]=(X[i].real() * cos(i*PI/(count*2)) + X[i].imag() * sin(i*PI/(count*2))) * dTemp;
}

// 释放内存
delete X;
}/*************************************************************************
 *
 * 函数名称：
 *   DIBIFFT()
 *
 * 参数:
 *   LPSTR lpDIBBits    - 指向源DIB图像指针
 *   long  lWidth       - 源图像宽度（象素数）
 *   long  lHeight      - 源图像高度（象素数）
 *  complex<double>*FD - 频域图象
 *  complex<double>*TD - 时域图象
 * 返回值:
 *   BOOL               - 成功返回TRUE，否则返回FALSE。
 *
 * 说明:
 *   该函数用来对经过FFT变换的图象进行反变换
 *
 ************************************************************************/
bool DIBIFFT(BYTE * lpDIBBits,long lWidth,long lHeight,
complex<double>*FD,complex<double>*TD)
{
// 指向源图像的指针
unsigned char* lpSrc;

// 中间变量
double dTemp;

// 循环变量
long i;
long j;

// 进行付立叶变换的宽度和高度（2的整数次方）
long w;
long h; double dMin,dMax;
long lCount=0;

int wp;
int hp;

// 图像每行的字节数
long lLineBytes;

// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;

// 计算进行付立叶变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
} // 计算图像每行的字节数
lLineBytes = WIDTHBYTES(w * 8);

FreqChange(FD,w,h);
for(i = 0; i < h; i++)
{
// 对y方向进行快速反付立叶变换
IFFT(&FD[w * i], &TD[w * i], wp);
}
// 保存变换结果
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
FD[ j*h+i ] = TD[ i*w+j ];
}
} for(i = 0; i < w; i++)
{
// 对x方向进行快速反付立叶变换
IFFT(&FD[i * h], &TD[i * h], hp);

}
/*******************倒序*********************/
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
FD[(h-1-i)*w+j] =TD[ i+h*j ];
}
}

for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
TD[ i*w+j ] = FD[ i*w+j ];
}
}
/******************倒序*************************/ dMin = dMax = TD[0].real();
//应该显示实部
for( i=0;i<w*h;i++ )
{
if( dMin>TD[i].real() )
dMin = TD[i].real(); if( dMax<TD[i].real() )
dMax = TD[i].real();

} lCount = long (dMax - dMin);

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{

//dTemp = (TD[(h-1-i)*w+j].real()-dMin)/lCount*255.0;
//dTemp = TD[(h-1-i)*w+j].real();
// 判断是否超过255
dTemp = abs(TD[(h-1-i)*w+j].real());
if (dTemp > 255)
{
// 对于超过的，直接设置为255
dTemp = 255;
}


// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)个象素的指针
// 此处不直接取i和j，是为了将变换后的原点移到中心
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * 
( h-1-i)+j; // (h - 1 - (i<h/2 ? i+h/2 : i-h/2)) 
// + (j<w/2 ? j+w/2 : j-w/2);

// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
} // 返回
return TRUE;
}

Answer 2

/*************************************************************************
*
* 函数名称：
*   DIBFFT()
*
* 参数:
*   LPSTR lpDIBBits    - 指向源DIB图像指针
*   long  lWidth       - 源图像宽度（象素数）
*   long  lHeight      - 源图像高度（象素数）
*  complex<double>*FD - 频域图象
*  complex<double>*TD - 时域图象
* 返回值:
*   BOOL               - 成功返回TRUE，否则返回FALSE。
*
* 说明:
*   该函数用来对时域图象进行FFT变换
*
************************************************************************/
bool DIBFFT(BYTE * lpDIBBits,long lWidth,long lHeight,
complex<double>*TD,complex<double>*FD)
{
// 指向源图像的指针
unsigned char* lpSrc;

// 中间变量
double dTemp;

// 循环变量
long i;
long j;

// 进行付立叶变换的宽度和高度（2的整数次方）
long w;
long h;

int wp;
int hp;

// 图像每行的字节数
long lLineBytes;



// 赋初值
w = 1;
h = 1;
wp = 0;
hp = 0;

// 计算进行付立叶变换的宽度和高度（2的整数次方）
while(w * 2 <= lWidth)
{
w *= 2;
wp++;
}

while(h * 2 <= lHeight)
{
h *= 2;
hp++;
} // 计算图像每行的字节数
lLineBytes = WIDTHBYTES(w * 8);

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 指向DIB第i行，第j个象素的指针
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;

// 给时域赋值
TD[ i*w+j ] = complex<double>(*(lpSrc), 0);
}
}

for(i = 0; i < h; i++)
{
// 对y方向进行快速付立叶变换
FFT(&TD[w * i], &FD[w * i], wp);
}

// 保存变换结果
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
TD[i + h * j] = FD[ i*w+j ];
}
}

for(i = 0; i < w; i++)
{
// 对x方向进行快速付立叶变换
FFT(&TD[i * h], &FD[i * h], hp);
}/*******************倒序**********************/
for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
TD[i*w+j] = FD[ i+h*j ];
}
}

for(i = 0; i < h; i++)
{
for(j = 0; j < w; j++)
{
FD[ i*w+j ] = TD[ i*w+j ];
}
}
/*******************倒序*************************/
FreqChange(FD,w,h);

// 行
for(i = 0; i < h; i++)
{
// 列
for(j = 0; j < w; j++)
{
// 计算频谱
dTemp = abs( FD[ (h-1-i)*w+j ] ) / 100;

// 判断是否超过255
if (dTemp > 255)
{
// 对于超过的，直接设置为255
dTemp = 255;
}

// 指向DIB第(i<h/2 ? i+h/2 : i-h/2)行，第(j<w/2 ? j+w/2 : j-w/2)个象素的指针
// 此处不直接取i和j，是为了将变换后的原点移到中心
//lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (lHeight - 1 - i) + j;
lpSrc = (unsigned char*)lpDIBBits + lLineBytes * 
( h-1-i)+j; // (h - 1 - (i<h/2 ? i+h/2 : i-h/2)) 
// + (j<w/2 ? j+w/2 : j-w/2);

// 更新源图像
* (lpSrc) = (BYTE)(dTemp);
}
} // 返回
return TRUE;}
这是我根据VC图象处理这本书改的，
包括了2维反FFT，
没错的！

Answer 3

matlab里有fft的函数可直接使用