二维傅里叶变换
上接文章
第一节 二维傅里叶变换对
在之前的章节所讨论的都是一维离散信号的傅里叶变换,如果将一维拓展到二维上,那么冲击采样函数应该满足如下描述:
同时,对二维信号的采样可以写成式7.2的形式
\int_{x_{0} = - \infty}^{+ \infty}{\int_{y_{0} = - \infty}^{+ \infty}{f(x,y)\delta\left( x - x_{0},y - y_{0} \right)\text{dxdy}}} = f\left( x_{0},y_{0} \right)\ (7.2)
如果采样的信号是离散的,那么,对二维信号的式子就应该由积分变为累加,对其采样如7.3所示:
\sum_{x_{0} = - \infty}^{+ \infty}{\sum_{y_{0} = - \infty}^{+ \infty}{f(x,y)\delta\left( x - x_{0},y - y_{0} \right)}} = f(x_{0},y_{0})
那么依据上述推论,就有理由将傅里叶变换的二维形式写作如7.3式:
F\left( u,v \right) = \sum_{x = 0}^{M - 1}{\sum_{y = 0}^{N - 1}{f(x,y)}}e^{\frac{- 2\text{πj}}{M}ux + \frac{- 2\text{πj}}{N}\text{vy}}\ (7.3)
f\left( x,y \right) = \frac{1}{\text{MN}}\sum_{u = 0}^{M - 1}{\sum_{v = 0}^{N - 1}{F(u,v)}}e^{\frac{2\text{πj}}{M}\text{ux} + \frac{2\text{πj}}{N}\text{vy}}\ (7.4)
7.3,7.4即为二维傅里叶变换对,其中离散信号F(u,v)是一个长为M,宽为N的二维矩阵,7.3位其正变换式,7.4位其逆变换式子,可以看到二维傅里叶变换实际上是一维傅里叶变换在横纵方向的卷积。
第二节 用C语言实现二维傅里叶变换对
从第一节的推论得知二维傅里叶变换实际上是一维傅里叶变换在横纵方向的卷积。因此，二维傅里叶变换可以在第六章的代码的基础上进行进一步拓展。
实现代码如下所示：
void FFT2(complex x[],complex X[],int N)
//对矩阵每一行进行快速傅里叶变换
for (int i=0;i<N;i++)
FFT(&x[i*N],&X[i*N],N);
//矩阵转置
for (int cy=0;cy<N;cy++)
for (int cx=cy+1;cx<N;cx++)
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
//对矩阵的每一列进行快速傅里叶变换
for(int i=0;i<N;i++)
FFT(&X[i*N],&X[i*N],N);
//再次进行矩阵转置,得到最终的变换结果
for (int cy=0;cy<N;cy++)
for (int cx=cy+1;cx<N;cx++)
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
其中,二维FFT的函数原型及实现如上所示,x[]表示输入的离散二维信号,X[]表示输出的信号,N表示该矩阵的边长,使用该函数的限制条件是,输入的二维信号矩阵必须是长宽相等,同时,边长满足基2数。该函数的实现原理是，首先输入一个二维矩阵（
在参数中表示为一维，即二维序列从左上角到右下角依次编号，如下所示
\begin{matrix} x_{0} & x_{1} & \ldots \ x_{n} & x_{n + 1} & \ldots \ \ldots & \ldots & \ldots \ \end{matrix}
），对该矩阵的每一行进行一维傅里叶变换,然后将变换结果(也是一个二维矩阵)进行转置,再对转置后的矩阵每一行进行傅里叶变换,最后再对变换的矩阵进行一次矩阵转置,就得到了最终变换后的结果。
同样的，二维傅里叶逆变换代码如下：
void IFFT2(complex X[],complex x[],int N)
//矩阵转置
for (int cy=0;cy<N;cy++)
for (int cx=cy+1;cx<N;cx++)
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
//对矩阵的每一行进行一维傅里叶逆变换
for(int i=0;i<N;i++)
IFFT(&x[i*N],&X[i*N],N);
//再次进行矩阵转置
for (int cy=0;cy<N;cy++)
for (int cx=cy+1;cx<N;cx++)
complex _t=X[cy*N+cx];
X[cy*N+cx]=X[cx*N+cy];
X[cx*N+cy]=_t;
//对转置后的矩阵每一行再次进行傅里叶逆变换
for (int i=0;i<N;i++)
IFFT(&X[i*N],&X[i*N],N);
到此,二维傅里叶变换对的代码也就编写完成了。
第三节 验证FFT2/IFFT2代码
为了验证上一节傅里叶变换对的代码准确性,编写测试代码对其进行测试:
int main()
//构造一个长度为8的复信号样本,与一个4x4的二维样本序列
complex samples[8],_out[8],samples_2d[4*4],_out2d[4*4];
int i,j;
//初始化这个复信号样本序列为1,2,3,4,5,6,7,8
for (i=0;i<8;i++)
samples[i].re=i;
samples[i].im=0;
//对该信号进行DFT运算
printf("DFT:\n");
//输出运算后的结果
DFT(samples,_out,8);
for (i=0;i<8;i++)
printf("(%f,%f)\n",_out[i].re,_out[i].im);
//对该信号进行IDFT逆变换
printf("IDFT:\n");
IDFT(_out,samples,8);
//输出结果,应该与原信号相一致
for (i=0;i<8;i++)
printf("(%f,%f)\n",samples[i].re,samples[i].im);
//重新初始化原信号(因为经过逆变换存在精度误差)
for (i=0;i<8;i++)
samples[i].re=i;
samples[i].im=0;
//进行FFT变换,应该与DFT变换的结果一致
printf("FFT:\n");
FFT(samples,_out,8);
//输出结果
for (i=0;i<8;i++)
printf("(%f,%f)\n",_out[i].re,_out[i].im);
//进行IFFT逆变换,应能够还原出原信号
printf("IFFT:\n");
IFFT(_out,samples,8);
//输出结果
for (i=0;i<8;i++)
printf("(%f,%f)\n",samples[i].re,samples[i].im);
//初始化二维样本,从左上角到右下角依次为0~16
for (i=0;i<16;i++)
samples_2d[i].re=i;
samples_2d[i].im=0;
printf("FFT2:\n");
//进行二维快速傅里叶变换
FFT2(samples_2d,_out2d,4);
//输出变换结果
for (i=0;i<4;i++)
for (j=0;j<4;j++)
printf("(%f,%f)",_out2d[i*4+j].re,_out2d[i*4+j].im);
printf("\n");
printf("\n");
//二维快速傅里叶逆变换
printf("IFFT2:\n");
//正变换的结果作为输入
IFFT2(_out2d,samples_2d,4);
//输出结果
for (i=0;i<4;i++)
for (j=0;j<4;j++)
printf("(%f,%f)",samples_2d[i*4+j].re,samples_2d[i*4+j].im);