《信号处理原理》实验指导

《信号处理原理》实验指导
【实验数据】
实验数据为一图像数据文件，文件格式是纯文本格式。

文件正文的第一行的值表示矩阵的大小，即N值。

后面的N行是点阵图像，每行有N个数据。

N最大为256。

在图像点阵中，‘.’代表0（即没有点），‘o’代表1（即有点）。

【实验要求】
(1)对输入图象文件内容进行2D－FFT变换，再对所得频谱数据进行
2D－IFFT，将结果重新转换成字符文件保存起来。

(2)对输入图象文件内容进行2D－FFT变换，再将频谱的大小压缩为
（N/2）*（N/2），然后对所得频谱数据进行2D－IFFT，将结果重
新转换成字符，保存成（N/2）*（N/2）的汉字图像。

(3)对输入图象文件内容进行2D－FFT变换，再将频谱的大小压缩为
（N/2）*（N/2），然后对频谱补入一些零，再对频谱进行2D－IFFT，
将结果重新转换成字符，保存成（N/2）*（N/2）的汉字图像。

【实验原理】
一幅二维数字图像可以用矩阵[g(m,n)]来表示，g(m,n)是图像在坐标(m,n)处的灰度级（或彩色RGB值）。

也可以把g(m,n)视为一个二元函
数，它的自变量为m 和n ，则可以用它来表示数字图像在平面上的亮度分布。

矩阵可以写成下面的形式：
[](0,0)
(0,1)(0,1)(1,0)(1,1)(1,1)[](,)(1,0)(1,1)
(1,1)g g g N g g g N g g m n g M g M g M N -⎡⎤⎢⎥-⎢⎥
==⎢⎥⎢
⎥
----⎣
⎦
在上面的基础上，我们可以定义下面的二维DFT ： 定义1：二维矩阵向量[g(m,n)]的2D-DFT
[]101
0(,)))(,,(N M m M m n p q N n G p q g m n W F g m W n --==⎡⎤⎢⎥⎣⎦=⎡⎤⎥⎣⎦
=⎢∑∑01p M ≤≤-, 01q N ≤≤-
从上面的定义我们可以看出：2D-DFT 可以用两次1D-DFT 来实现。

一次是对g(m,n)的各行（即m 相同而n 不同）进行1D-DFT---如上式中的红色部分所示；再对变换后的结果，再按列进行1D-DFT---如上式中的蓝色部分所示。

上面这个公式是我们实现2D-DFT 的算法基础。

定义2：二维矩阵谱向量[G(p,q)]的2D-IDFT
[]1
1
1
1(,1
(,)(,))N nq N q M mp
M p G p q g m n F G p M
W q N
W ----=-===
∑∑
01m M ≤≤-,01n N ≤≤-
根据相同的想法，上面的2D-IDFT 公式也可以用两次1D-IDFT 来实现。

【2D-FFT 算法实现】
二维FFT 相当于对行和列分别进行一维FFT 运算。

具体的实现办法如下： 先对各行逐一进行一维FFT ，然后再对变换后的新矩阵的各列逐一进行一维FFT 。

相应的伪代码如下所示：
for (int i=0; i<M; i++)
FFT_1D(ROW[i]，N); for (int j=0; j<N; j++)
FFT_1D(COL[j]，M);
其中，ROW[i]表示矩阵的第i 行。

注意这只是一个简单的记法，并不能完全照抄。

还需要通过一些语句来生成各行的数据。

同理，COL[i]是对矩阵的第i 列的一种简单表示方法。

所以，关键是一维FFT 算法的实现。

下面讨论一维FFT 的算法原理。

【1D-FFT 的算法实现】
设序列h(n)长度为N ，将其按下标的奇偶性分成两组，即h e 和h o 序列，它们的长度都是N/2。

这样，可以将h(n)的FFT 计算公式改写如下 ：
n
n
()()()nk nk N
N
H k h n W
h n W
=
+
∑∑偶数奇数
()
()
2
2
1
1220
(2)(21)N N rk rk k N
N
N
r r h r W
W
h r W --===++∑∑………………… (A)
由于
2/2N N W W =
所以，(A)式可以改写成下面的形式：
22
2
2
1
10
()()N N N N rk k rk
N r r e o W W h r h r W --===+∑∑ 按照FFT 的定义，上面的式子实际上是：
()()k
e N o H k W H k =+
其中，k 的取值范围是 0~N-1。

我们注意到He(k)和Ho(k)是N/2点的DFT ，其周期是N/2。

因此，H(k)DFT 的前N/2点和后N/2点都可以用He(k)和Ho(k)来表示
22()()()()
N e e N o o H k H k H k H k +=+=
而且
22
()N N k k k
N
N N N W
W W W +=⋅=-
于是，N 点H(k)用N/2点的He(k)和Ho(k)来计算的公式为：
()()()k e N o H k H k W H k =+
2()2
2
2()()()N k N
N N e N
o H k H k W
H k ++=+++
()()2k e N o N H k H k W H k ⎛⎫
+=- ⎪⎝⎭
20,1,
,1N k =-
【1D-FFT 的算法流程】
根据上面推导的公式，我们可以用递归程序来实现1D-FFT 。

下面是一个示例：
void FFT_1D (Comp in[ ], Comp out[ ], int N)
{
Comp he[256], ho[256];
Comp He[256], Ho[256];
// 如果DFT点数为1，则根据DFT公式，可以直接返回原值
// 这是递归的退出条件
if (N==1)
{
out[0] = in[0];
}
else { // 如果N不是1，则…
// 按下标将数据分成两组
for (int i=0; i<N/2; i++)
{
he[i] = in[2*i]; // 偶数下标分组
ho[i] = in[2*i+1]; // 奇数下标分组
}
// 计算偶数下标部分的N/2点DFT
FFT_1D(he, He, N/2);
// 计算奇数下标部分的N/2点DFT
FFT_1D(ho, Ho, N/2);
// 计算当前所有数据的N点DFT
for (int k=0; k<N/2; k++)
{ // 下面计算式中的W[k]是旋转因子，可以事先计算好。

// 当前数据的前一半，N/2点
out[k] = He[k] + W[k]*Ho[k];
// 当前数据的后一半，N/2点
out[N/2+k] = He[k] –W[k]*Ho[k];
}
}
}
【图象压缩与放大原理】
要尽量少失真地对图象进行缩放处理，则需从频域着手进行。

如果从频谱图中去掉一些频率的分量，然后再作IFFT，得到时域数据，即可实现对图象的缩小。

而如果在频谱中加入一些零值，再作IFFT，得到时域数据，则可
实现对图象的放大。

根据一维情形下的奈奎斯圣区间和折叠频率的概念，我们知道，二维频谱的高频区集中在二维谱矩阵（谱向量的矩阵形式）的中心区域；而低频区在4个角上。

如果考虑到二维频谱图中高频分量与低频分量的分布情况，则对图象频谱图的处理，不论是去掉部分频率数据，还是在频谱数据中补零，都应关于图象频谱的中心对称。

要思考和解决的问题是：为保证图象缩小后失真最小，应去掉原图中的高频分量，还是去掉低频分量。

【注意事项】
(1)输入的图像数据是以文本形式存放在文件中的，在读入文件内容后要
将不同的字符转化为数值（如·转为0，o转为10），才能进行FFT
变换。

(2)为了能显示还原后的图像，还需要将IFFT变换后的数据（复数值）
转换成不同的字符，并以文本形式入文件中。

(3)复数转成字符的原则和方法：IFFT变换后的数据是复数，根据复数
模的大小，将它们分别转成·和o字符。

这就要设定一个控制阈值，
模大于阈值的复数对应o字符，模小于阈值的复数对应字符。

这个控
制阈值可能通过实验来调整和确定。

(4)程序中涉及到了复数运算（加减乘除），可以自己编程实现，也可以
通过包含下面的代码来解决：
#include <complex>
using namespace std;
// 将complex<float>简记为Comp
typedef complex<float> Comp;
这样，就可以在程序中直接使用“Comp”这种类型了。

用法如下：Comp w1(3, 4);
Comp w2;
W2 = Comp(3, 4);
【主程序代码流程示例】
(1)接受用户对阈值的设定
(2)打开文件，读入字符
(3)根据（1）中设定的阈值，将字符转换成数值保存到数组
(4)调用2D-FFT，计算图象的频谱
(5)对频谱数据进行相应处理（如去掉一部分频谱值，或添加一些零）
(6)对处理后频谱数据进行2D-IFFT，转换成时域数值
(7)根据阈值，将还原的时域数据转换成字符
(8)将字符保存到文件，以备查看。