基于FPGA的数字图像处理-线性滤波器【3.6】

式中，D0 为滤波器的截止频率，n为滤波器的阶数，W为滤波器的 带宽。 以4阶巴特沃斯滤波器为例，C语言示例代码如下： //计算巴特沃斯高通滤波器函数（4阶）

double Cal_Butterworth_HPF(int u,int v,double sigma)
{
double temp = 0;
temp=1+((sigma*sigma)/(u*u+v*v))*
((sigma*sigma)/(u*u+v*v));
temp = 1/temp;
return temp;
}
//计算巴特沃斯低通滤波器函数（4阶）
double Cal_Butterworth_LPF(int u,int v,double sigma)
{
double temp = 0;
temp = 1 + ((u*u + v*v)/(sigma*sigma))*((u*u +
v*v)/(sigma*sigma));
temp = 1 / temp;
return temp;
}
//计算二维巴特沃斯传递函数：带阻滤波器（4阶）
double Cal_Butterworth_BSF(int u,int v,double
sigma,double width)
{double temp=((u*u+v*v)-(sigma*sigma))*((u*u+v*v)-
(sigma*sigma));
double temp1 = (u*u + v*v)*width*width;
double temp2=(temp1/temp)*(temp1/temp)*(temp1/temp)*
(temp1/temp);
temp2 = 1 / (1 + temp2);
return temp2;
}
//计算二维巴特沃斯传递函数：带通滤波器（4阶）
double Cal_Butterworth_BPF(int u,int v,double
sigma,double width)
{
double temp=((u*u+v*v)-(sigma*sigma))*((u*u+v*v)-
(sigma*sigma));
double temp1 = (u*u + v*v)*width*width;
double temp2=(temp1/temp)*(temp1/temp)*(temp1/temp)*
(temp1/temp);
temp2 = 1 - 1 / (1 + temp2);
return temp2;
}

下面的C++代码展示了如何利用FFT来实现一副图像的低通滤波功能

Transform_FFT *p_fft_orig = new Transform_FFT; //定义原图
FFT变换对象
BYTE* p_FFT_orig = new BYTE[m_dwHeight*m_dwWidth];//原图
FFT结果p_fft_orig-
>ImgFFT2D(m_pTemp,m_dwHeight,m_dwHeight,p_FFT_orig);
//原图进行快速傅里叶变换
Transform_FFT *p_FFT_filt = new Transform_FFT;/*滤波后的
FFT结果*/
p_FFT_filt->m_pFFTBuf = new ComplexNumber[m_dwWidth*
m_dwHeight];
for (i = 0; i＜m_dwHeight; i++)
for (j = 0; j＜m_dwWidth; j++)
{
if(i＜m_dwHeight/2)i0=i+m_dwHeight/2;
else i0=i-m_dwHeight/2;
if(j＜m_dwWidth/2)j0=j+m_dwWidth/2;
else j0=j-m_dwWidth/2; //将频谱搬移到图像中心
/*计算当前频谱值*/
a = p_fft_orig->m_pFFTBuf[i0*m_dwWidth + j0].real;
b = p_fft_orig->m_pFFTBuf[i0*m_dwWidth + j0].imag;
//计算巴特沃斯低通滤波器增益
Gaus_Value = Cal_Butterworth_LPF(i - m_dwWidth / 2,j -
m_dwWidth /2,70);
p_FFT_filt->m_pFFTBuf[i0*m_dwWidth + j0].real =
a*Gaus_Value;
/*频域直接相乘*/
p_FFT_filt->m_pFFTBuf[i0*m_dwWidth + j0].imag =
b*Gaus_Value;
}
BYTE *p_ifft_fusion = new BYTE[m_dwHeight*m_dwWidth];//对滤波后结果极性fft反变换
p_FFT_filt-
>ImgIFFT2D(p_ifft_fusion,m_dwWidth,m_dwHeight);

运算结果如图7-6（b）所示，图7-6（a）为原始图像，图7- 6（b）为滤波后的图像。由图7-6可以看出，频率域的低频滤波和空域 中的低频滤波可以达到相同的效果，不同的是时域需要卷积，而频域 是直接相乘。

图7-7是对一个被固定频带干扰的图像经过4阶巴特沃斯带阻滤波 器滤波后的效果图。图7-7（a）为原图，图7-7（b）为处理过后的图 像。

由处理后的图可见，高频信息已经被滤除。不妨对比一下两幅图 的FFT频谱，如图7-8所示，图（a）为被干扰的频谱，图（b）为处理 后的频谱。

所设计的带阻滤波器有效滤除了干扰频段，但同时滤除了原图中 的一部分高频信息，因此相对于原图稍微模糊。同样，对干扰过的原图在时域进行均值滤波，结果如图7-9（b） 所示。

        可以看到，处理过后的图像还是有不少的高频纹理。如果加大滤 波窗口尺寸，高频纹理就可以在一定程度上有所缓解，但是图像丢失 的细节会更多。 7.2 基于FPGA的均值滤波 正如在上一节所讲到的，线性滤波实现简单，十分适合用FPGA来 实现。一般情况下，FPGA在前端捕获到视频数据之后首先需对图像数 据做一个简单的预处理，然后根据噪声的来源，针对椒盐噪声进行中 值滤波处理，针对高斯噪声进行高斯滤波处理，均值滤波在预处理中 也十分常见。同时，边缘提取及梯度计算也是许多复杂处理算法的基 础。下文将详细介绍均值滤波器以及Sobel算子的FPGA实现。

7.2.1 整体设计与模块划分

        再把均值滤波的数学表达式列出如下：

由上述公式列出求图像均值的步骤： （1）获得当前窗口所有像素。 （2）计算当前窗口所有像素之和。 （3）将第（2）步结果除以当前窗口数据总数。 （4）滑动窗口到下一个窗口，直到遍历完整幅图像。 滤波采用滑动窗口方法来实现整幅图像的遍历，因此，采用流水 线结构来设计是再也合适不过的了。对于流水线结构来说，每个像素 的运算方法是一致的，所需考虑的只是边界像素的处理问题。前面已 经详细介绍了如何对图像进行行列寻址及如何进行行列对齐，因此， 本模块的设计重点在计算窗口像素和和除法操作。 以5×5的均值滤波窗口为例，顶层设计如图7-10所示。

首先来看二维窗口求和模块。 一般情况下，任何二维的计算步骤都可以化为一维的操作。我们 在前面的对行列对齐的介绍中也提到，对于图像来说，就是一个二维数组。由于行方向的数据流是连续的，因此在流水线操作中，常常会 首先进行行方向的操作。 假定现在已经完成了第一行的求和操作，接下来需要“等”下行 的求和完成。如何进行等待？在FPGA中，等待的实现方法就是进行缓 存。二维操作转换为一维操作后的结构如图7-11所示。

接下来的问题是，如何进行一维向量求和操作？对于1×5的向量 求和而言，当前数据需要“等到”下4个数据到来之后才能得到连续5 个数据，并执行加法操作。可以预期的是，还是需要把前几个数据单 独缓存起来，一个指定位宽的寄存器即可满足要求。同步5个连续的输 入数据如图7-12所示。

最后的问题是求取窗口的均值，需要将上述计算出来的和除以一 个归一化系数，也就是整个窗口的像素数目。正如前面所讲到的，在 FPGA里面对于一方确定的除法操作，一般情况下不直接进行除法运 算，而是通过近似的乘加方法来实现等效转换。在这里，对于固定的窗口，除法的分母是固定的。因此，完全可 以用此方法来实现等效近似。具体的转换方法将在下一节进行介绍。 7.2.2 子模块设计 按照上一节的整体设计，需要设计以下几个子模块：

（1）一维求和模块，这里记为Sum_1D。

（2）二维求和模块这里记为Sum_2D。

（3）除法转换模块，此模块比较简单，一般情况下不进行模块封 装。

（4）行缓存电路实现行列间像素对齐。 整个顶层模块调用Sum_2D模块和除法转换电路来实现求取均值， 记为Mean_2D。 1.一维求和模块设计（Sum_1D） 用FPGA来求和是再简单不过的操作了，求和操作也是FPGA所擅长 的事情之一。所要注意的只是求和结果不要溢出。一般情况下，2个位 宽为DW的数据相加，至少得用一个DW+1位宽的数据来存放。 由于是求取连续数据流的和，正如前面所讲到的，最简单的办法 是将数据连续打几拍，对齐后进行求和。假定窗口尺寸为5，则求和电 路可以根据图7-13所示进行设计（读者可以思考一下为什么不直接将5 个数相加）。

图7-13所示的设计思路非常简单，将输入数据流连续打4拍，加上 当前数据组成连续5拍数据，经过3个时钟的两两相加运算，即可得到 连续5个数据的和。 由图7-13可知，此电路的资源消耗为4个加法器、7个寄存器，运 算开销为3个时钟。 当然上面的电路确实可以实现预定的功能，然而本书中采用另外 一种方法。这种方法就是利用增量更新的方式来实现窗口横向求和， 这种求和方式在大尺寸的窗口计算中十分有用。 在连续两个像素求和的过程中，仅仅有头尾的两个像素不同。假 定当前计算地址为n+1，计算结果为Sum(n+1)，上一个地址为n，计算 结果为Sum(n)，输入数据流为X(i)，设定当前计算窗口尺寸为7，则有

也就是针对每一个窗口并不需要重新计算所有窗口内的像素和， 可以通过前一个中心点的像素和再通过加法将新增点和舍弃点之间的 差计算进去就可以获得新窗口内像素和。 具体到FPGA实现方面，同样需要把数据连续打几拍，同时计算首 个数据与最后一个数据的差。当前求和结果为上一个求和结果与计算 之差的和。同样对于窗口尺寸为5的行方向求和操作，设计电路如图7- 14所示。 由图7-14可知，此电路只需1个加法器和1个减法器。可以预见的 是，无论窗口尺寸多大，所需的加法器和减法器也都是1个。因此，在 窗口尺寸比较大的情况下，可得到比第一个设计电路更优的资源消耗 的目的。不仅如此，求和电路的计算开销仅为1个时钟。

2.二维求和模块设计（Sum_2D） 目前我们已经实现了窗口内一维行方向上的求和工作，现在要得 到整个窗口内的像素之和，还必须将每一行的计算结果再叠加起来。 那么每一行的计算结果是否也可以采取上面的增量更新的方法进行计 算？答案显然是否定的，这是由于纵向的数据流不是流水线式的。这 时，就必须要采用第一种方法(见图7-13)所采用的求和方式。 同样，在进行列方向上的求和时，需要进行行缓存，并将一维行 方向的求和结果打入行缓存，行缓存的个数为窗口尺寸减去1。 就窗口尺寸5×5的情况而言，二维求和模块的电路设计如图7-15 所示。

3.除法电路设计 在前面也已经提到，对于分母固定的除法操作，可以通过泰勒展 开或是移位转换等方式转换为FPGA所擅长的移位、加法与乘法操作。 还是首先假定窗口尺寸为5，在求取窗口像素和之后需要除以25来 求得均值。假定求得结果为Sum，计算后的均值为Average，则有 

实际上，在计算的过程中，可以通过控制最终相加结果的移位位 数来保留小数位，从而提高计算精度。例如，若想要保留3位小数位， 则上面的公式就变为

此时将结果右移7位即可达到我们的要求。这样，便将以此除法操 作转换为9次移位操作和7次加法操作。实际上，对于FPGA来说，移位 操作是0开销（不消耗时钟周期）的。同样，对于其他尺寸的处理窗口，也可以采用类似的转换方法。 以5×5的窗口为例，将除法电路加上之后得到的求均值电路如图7-16 所示。

7.2.3 Verilog代码设计

        在代码设计中，其中一个需要考虑的问题是可重用性。在模块设 计中，建议将一些参量参数化，例如数据的位宽、图像的宽度和高 度，以及处理窗口的宽度等。对于求窗口均值操作来说，所遇见的问 题是，对于不同的窗口尺寸，除法电路是完全不一致的，我们所采用 的思路是用verilog的generate语句来实现条件编译的。 以下将详细介绍各个模块的代码设计过程。 1.一维求和模块设计（Sum_1D） 由图7-14可知，一维求和电路的设计十分简单，仅需若干个寄存 器和一个加法器和1个减法器即可满足要求。设计的重点在于前端数据 寄存器的设计。 不妨将所采用的处理窗口尺寸（KernelSize）记为KSZ，数据位宽 记为DW，也将这两个参数作为可配置形参。模块定义如下：

module sum_1d(clk, //同步时钟
din, //输入数据流
din_valid, //输入数据有效
dout_valid, //输出数据有效
dout //输出数据流
);
parameter DW = 14; //数据位宽参数
parameter KSZ = 3; //求和窗口参数
//定义KSZ+1个输入寄存器
reg [DW-1:0] reg_din[0:KSZ];
//定义上一个求和寄存器
reg [2*DW-1:0]sum;
//定义中间信号
wire [2*DW-1:0]sub_out;
//定义减法器输出信号
wire [2*DW-1:0]diff;
//连续缓存KSZ拍信号 同时缓存输入有效信号
always @(posedge clk)
begin
din_valid_r ＜= #1 ({din_valid_r[KSZ -
1:0],din_valid});
reg_din[0]＜= #1 din;
for (j = 1; j ＜= KSZ; j = j + 1)
reg_din[j]＜= #1 reg_din[j - 1];
end
//做减法计算差值assign sub_out=
((din_valid_r[0]==1'b1&din_valid_r[KSZ]==1'b1)) ?
({{DW{1'b0}},reg_din[KSZ]}) :
({2*DW{1'b0}});
assign diff = ({{DW{1'b0}},reg_din[0]}) - sub_out;
//计算最后的求和结果
always @(posedge clk)
begin
if (din_valid == 1'b1 & ((～(din_valid_r[0]))) == 1'b1)
sum ＜= #1 {2*DW-1+1{1'b0}};
else if ((din_valid_r[0]) == 1'b1)
sum ＜= #1 sum + diff;
end
//输出信号
assign dout_valid = din_valid_r[1];
assign dout = sum;