离散小波变换的FPGA-Verilog实现

to 51研究不顺的假期

[toc]

在上一次的博客中我们介绍了短时傅里叶变换,小波函数,尺度函数,离散小波变换.详情请看:

小波变换(一)

这次我们把把整个变换说完,顺便介绍一下他的硬件实现(一维)

尺度函数族

众所周知,小波变换的双正交基就来自与小波函数和尺度函数,而他们通过scale和平移来得到的小波函数族和尺度函数族正正表示了不同小波(尺度)函数的分辨率,下面来用公式回顾一下:

小波函数定义

小波函数族$\psi_{j,k}(t)$定义为:
$$\psi_{j,k}(t) = 2^{j/2}\psi(2^jt-k);\quad j,k\in Z$$

尺度函数定义

尺度函数族$\varphi_{j,k}(t)$定义为:
$$\varphi_{j,k}(t) = 2^{j/2}\varphi(2^jt-k);\quad j,k\in Z$$

尺度函数分辨率沿拓

事实上,尺度函数族中的尺度 j 决定着不同的时间分辨率

而实际上(证明请看参考资料),高时间分辨率的尺度函数必然可以代表低时间分辨率的尺度函数,如上图的三角尺度函数可以表示为:

多分辨分析(MRA)方程

尺度函数的MRA方程

有了上面的铺垫,这里就比较简单了,因为我们知道,低分辨率信号可以由高分辨率信号线性表达,所以我们可以看出,高分辨率信号所张成的空间必然包含低分辨率张成的空间:

图中越外围的代表越高分辨率的尺度函数,而$ L^2 $为平方可积空间,即:
$$ \int_{t\in R}|x(t)|^2 dt < \infty $$
学过信号与系统的应该都知道积分式的定义吧…

所以在这个基础上我们引出尺度函数的多分辨分析方程(用以映射用高分辨率信号来表示低分辨率信号之间的关系)
$$
\varphi(t)=\sum_{k} h_{0}[k] \cdot \sqrt{2} \varphi(2 t-k)
$$
其中$h_0[k]$是尺度函数系数,也可以说成尺度滤波器的单位脉冲响应,到后面这个是重点.

小波函数的MRA方程

讲小波函数的MRA方程之前,我们不妨来回顾一下IDWT的定义:
$$x(t) = \sum_nc_0[k]\varphi_{0,k}(t) + \sum_nd_0[k]\psi_{0,k}(t) + \sum_nd_1[k]\psi_{1,k}(t) +\cdots $$

我们可以看到上面的IDWT中,我们给定了小波函数$varphi_{0,k}(t) $的尺度$j=0$ ,以此为基底,再用不同分辨率的尺度函数来描述整个x(t).

实际上,我们可以重构上式:
$$x(t)=\sum_{k} c_{j_{0}, k} \varphi_{j_{0}, k}(t)+\sum_{j=j_{0}}^{\infty} \sum_{k} d_{j, k} \psi_{j, k}(t)$$

假设,此时他的信号空间由小波函数族张成的$V_{j_0}$和尺度函数族张成的$$W_{j_0}\oplus W_{j_0 +1} \oplus W_{j_0 +2} \oplus W_{j_0+3} \oplus \cdots $$

则$L^2$空间可表示为:
$$ L^2 = V_{j_0} \oplus W_{j_0}\oplus W_{j_0 +1} \oplus W_{j_0 +2} \oplus W_{j_0+3} \oplus \cdots $$

而事实上,小波函数和尺度函数满足一下关系:
$$
\begin{array}{l}{V_{1}=V_{0} \oplus W_{0}} \ {x_{1}(t)=C_{0}(t)+D_{0}(t)}\end{array}
$$

所以上上式中可以将$ V_{j_0} $无限递归分解为小波信号张成的信号空间,即:
$$
\begin{array}{c}{L^{2}=\cdots \oplus W_{-2} \oplus W_{-1} \oplus W_{0} \oplus W_{1} \oplus W_{2} \oplus \cdots} \ {x(t)=\sum_{k} \sum_{j=-\infty}^{\infty} d_{j, k} \psi_{j, k}(t)}\end{array}
$$

所以从这个角度去考虑的话,我们就知道为什么$c_j [k]$表示的是信号的粗略(coarse)信息,而$ d_j [k] $ 表示的是信号的精细(detail)信息,当然,这只是其中一种理解,还有一种比较有意思的想法就是从差分的角度出发来考虑两种系数之间的意义,详见孙延奎教授的《小波分析及其应用》

结合尺度函数的MRA方程中的信号空间的图,我们可以看出,其实每一圈$V_j$所包含的就是低一维的尺度函数族张成的信号空间($V_{j-1}$)在加上低一维的小波函数张成的信号空间($W_{j-1}$),即:

自然而然,我们可以知道小波函数的MRA方程的递归意义是更为重要的,所以有:
$$
\psi(t)=\sum_{k} h_{1}[k] \cdot \sqrt{2} \varphi(2 t-k)
$$
其中$ h_1[k] $称为小波函数系数

由于篇幅关系,这里没办法进一步讨论信号空间映射,小波函数系数和尺度函数系数的关系还有对称性,看情况会再作一篇博客讨论.

离散小波变换的分解算法

有了MRA方程的概念,想必这里会比较简单,其实分解算法很简单,就是围绕着
$$
\begin{array}{l}{V_{1}=V_{0} \oplus W_{0}} \ {x_{1}(t)=C_{0}(t)+D_{0}(t)}\end{array}
$$
这个方程来进行分解,这里不给出证明过程(其实很简单,但是篇幅过长,详见参考资料),直接给出算法框图:
DWT的一级分解算法框图:

再将$C_j [k]$ 进行分解,可得DWT的二级分解算法框图:

所以多级的大家也会求了吧…

可以简单地看出,信号经过小波函数系数(尺度函数系数)之后还需要经过一个抽取的过程,这个的话自己看看小波函数和尺度函数的定义式就可以了

离散小波变换的重构算法

这里是相似的,而且由于懒的原因我没有做重构,所以也直接放图:

可以简单地看出,信号经过小波函数系数(尺度函数系数)之后还需要经过一个内插的过程,这个的话自己看看小波函数和尺度函数的定义式就可以了.

需要注意的是,可以先抽取再滤波,但是不能先滤波后内插,见上图,(常识)

小波变换(DWT)的FPGA实现

总所周知,这次我们要实现的算法框图是这个:

多相结构模型

本来我是直接按照算法流程实现了DWT,然后讲抽取和滤波对调了位置(为了提高系统的性能)

后来从网上仅有的资料中查看到别人做了多相分解,后面我想了想,如果不用多相结构的话,相当于原信号的先经过了一次抽取,也是极大地浪费了信号.
又有:
$$ H_0(z) = H_{0e} (z^2) +z^{-1} H_{0o}(z^2) $$
$$ H_1(z) = H_{1e} (z^2) +z^{-1} H_{1o}(z^2) $$
所以我们将信号和滤波器系数都进行奇偶分解,分别进行滤波,得到整体FPGA框图:

接下来我们简单地理一下过程

注意事项

激励和系数都是带符号位的,所有和信号流有关的reg和ip都要设置成signed
信号位宽是16,滤波器系数位宽是8
以db4小波为示例,其他小波只需要改改滤波器和matalb就可以了
仅实现一级分解,多级分解只要自己认真看博客就知道怎么做了

信号奇偶分解:

这个模块比较简单,但是想设计好需要一点小心思,思路就是做一个二分频时钟,上升沿将数据写入偶数部分,下降沿将数据写入奇数部分,还顺手设计了个三段式写法,确保建立时间和保持时间是信号保持时间的一半,但是都是FPGAer的常识了,这里不给出代码

matlab获取滤波器系数

代码很简单,先生成,再量化:

wn='db4';
[Ld,Hd,Lr,Hr] = wfilters(wn);
k=0:length(Ld)-1;
 subplot(221);stem(k,Ld);
 title('低通分解滤波器Ld');
 subplot(222);stem(k,Lr);
 title('低通重建滤波器Lr');
 subplot(223);stem(k,Hd);
 title('高通分解滤波器Hd');
 subplot(224);stem(k,Hd);
 title('高通重建滤波器Hr');

qua_ld = round(Ld*2^8);
qua_hd = round(Hd*2^8);
qua_lr = round(Lr*2^8);
qua_hr = round(Hr*2^8);
disp(qua_ld);
disp(qua_hd);
disp(qua_lr);
disp(qua_hr);

FPGA设计滤波器

这里的流程跟我上一篇博客,FPGA/Verilog 设计FIR滤波器是相似的,这里不多说,给出一个滤波器的源码:

module filter_hd_e(
	//clock and reset
	input CLK_50M,RST_N,
	//filter in out
	input  signed [15:0] data_in_hd_e,
	output signed [19:0] hd_out_e
);
localparam COEFF1,COEFF3,COEFF5,COEFF7;  //过长省略
wire signed [19:0] add_result;
reg signed [15:0] data_shift[3:0];
wire signed [23:0]  mul_data[3:0];	

 add_final add_inst(
	.clock       (CLK_50M),
	.data0x      (mul_data[0][23:8]),
	.data1x      (mul_data[1][23:8]),
	.data2x      (mul_data[2][23:8]),
	.data3x      (mul_data[3][23:8]),
	.result      (add_result)
);
always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)begin
        data_shift[0]  <= 0;
        data_shift[1]  <= 0;
        data_shift[2]  <= 0;
        data_shift[3]  <= 0;
    end
    else begin 
        data_shift[3]  <= data_shift[2];
        data_shift[2]  <= data_shift[1];
        data_shift[1]  <= data_shift[0]; 
        data_shift[0]  <= data_in_hd_e;
    end                  
end
mult	mult_inst_1 ( data_shift[0] *COEFF1 //下同
mult	mult_inst_2 (
mult	mult_inst_3 (
mult	mult_inst_4 ( //过长省略
assign hd_out_e = add_result;

endmodule

大家对应这框图可能会说我怎么没有把滤波器系数给翻转,这个问题的话,大家可以看看那个获取系数的函数描述,他本来就帮我们翻转了.

加法器就算了,过于简单

matlab产生激励

这个在上一篇博客中也有提及,但是这次我们不是直接产生mif文件,而是选择在仿真的时候读入数据,所以代码就是:

depth = 1024;
width = 16;
x = 0 : 2*pi/(depth-1) :2*pi;
y = sin(x)+sin(8*x);
y = (y/2) * 32768;%将信号放大32768倍
y(360:400) = 32767;    //为了小波变换的戏剧性效果设置
b = signed2unsigned(y,width);
fid = fopen('sinx.txt','wt'); %将信号写入一个.txt文件中
for num=0 : (depth-1) 
fprintf(fid,'%x\n',round(b(num+1)));
end
fclose(fid);

仿真结果

用tb读入数据:

integer i;   //数组坐标
reg signed [15:0] stimulus[1:data_num];  //数组形式存储读出的数据
initial 
begin
    RST_N = 1'b1;
	#60 RST_N = 1'b0;
	#60 RST_N = 1'b1; 
    $readmemh("sinx.txt", stimulus);  //将txt文件中的数据存储在数组中
    i = 0;
    repeat(data_num) begin   //重复读取数组中的数据
        i = i + 1;
        data_in = stimulus[i]; 
        #PERIOD;         //每个时钟读取一次
    end
	 $stop;
end