今日你学左米啊

0%

LMS自适应滤波的FPGA实现(二)

发表于 更新于 分类于 Valine：

LMS自适应滤波的FPGA实现(二)

闭环系统搞自闭著

[TOC]

算法流程

这个算法来源于上一篇博客LMS自适应滤波的FPGA实现(一)
所介绍的最后一个算法->Widrow-Hoff -LMS算法

回顾一下,最主要的公式:
$$ \vec{f}[n+1] = \vec{f}[n] + \mu e[n]\vec{x}[n] $$
其中:
$$ e[n] = d[n] - y[n] $$

所以我们可以大致得到算法流程:

  1. 初始化N阶向量 $$ \vec{f} = \vec{x} = \vec{0} $$
  2. 接收一组新的采样值{x[n],d[n]} ,并且将x[n]移入移位寄存器
  3. 计算FIR滤波器输出$$ y[n] = \vec{f^T}[n]\vec{x}[n] $$
  4. 更新误差函数$$ e[n] = d[n] -y[n] $$
  5. 更新滤波器系数$$ \vec{f}[n+1] = \vec{f}[n] + \mu e[n]\vec{x}[n] $$
  6. To Step 2

在FPGA中每一个Step单独用一个Always块描述就完事了.

应用场景

由于即便是一个最小的自适应滤波器系统也是一个闭环系统,所以考虑到调试的简便性,这里主要做的工作是干扰消除,是自适应滤波里面差不多最常用也是最简单的应用.同时也是一道国赛题.

在这个场景下,我们要实现的是,已知噪声信号和原始信号,需要从中还原出真实信号,所以系统框图是这样的:
又是一张自己做的图

FPGA实现4阶LMS自适应滤波

假设AD是12位,采样频率是50M的.
根据上面的讨论和假设,代码的大体框架其实已经给出了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
module fir_lms_one
#(parameter N1 = 12 , N2=25, L=4) 
	//N1 for Input data bit width and output coefficient bitwidth 
	//N2 for two times of N1 , +1 for 4 coef
	//L  for Filter length 
	//Delay for the pipeline of the system
(
	input clk,RST_N,					//system clk and reset
	input signed [N1-1:0] x_in,	//system_input
	input signed [N1-1:0] d_in,	//Reference input
	
	output signed [N1-1:0] f0_out, // 1 st coefficient
	output signed [N1-1:0] f1_out, // 2 nd coefficient
	output signed [N1-1:0] f2_out, // 3 rd coefficient
	output signed [N1-1:0] f3_out, // 4 th coefficient
	
	output signed [N2-1:0] y_out,	 //System_output
	output signed [N2-1:0] e_out	 //error sig
);

因为我的quartus装在ubuntu上面,大不了中文注释,大家大概看看…
接着我们只要一步一步跟着算法流程跑就完事了,这里有三个点:

  1. FIR滤波器的输出涉及到的乘法器计算
  2. 在更新误差函数和滤波器系数之间其实也涉及向量乘法
  3. 步进值的选择

下面就这三点给出部分源码(都很简单):

  1. FIR乘法器计算:
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    //the fir filter circuit
    reg signed [N1-1:0] f [0:3];
    reg signed [N2-1:0] p [0:3]; // fir cal temp array 
     always @(*) 
     begin : MulGen1    //for inst L multplier
       integer I;   		// loop variable 
       for (I=0; I<L; I=I+1) p[I] <= x[I] * f[I];
     end
     
    wire signed [N2-1:0] y;
    assign y = p[0]+p[1]+p[2]+p[3] ;		 //accumulation for fir

实际上,这就写完了…一定要注意这里的位宽对应,避免溢出

  1. 误差系数的更新和步进值
    1
    2
    3
    4
    5
    6
    7
    8
    9
    //the error cal circuit
    wire signed[N2-1:0] e;					  //to update the error 
    assign e 	= d_reg - (y>>>11);      //cause to the bitwidth  
    reg signed [N2-1:0] muex [0:3]; // error coef temp array 
    always @(*) 
     begin : MulGen2    //for inst L multplier
       integer I;    	// loop variable 
       for (I=0; I<L; I=I+1) muex[I] <= (x[I] * e)>>>4;
     end

这里的y因为是4个系数的乘积和,所以他的位宽和新的d[n]位宽是不一样的,需要截位处理

另外就是步进值的选择,实际上我见过很多人写的都是在算出e之后直接移位的,但实际上这种做法稍稍有点不妥.

在上一篇博客中,我们没有对这个不仅率进行讨论(虽然现在也不打算),但是,在理论上这一个值并不能太大,否则会导致梯度下降时在极值附近震荡,无法收敛.

除此之外,我们更需要考虑的一个问题是,在硬件实现时的字长效应.如果我们在算出误差后直接移位的话,很容易会导致误差直接截位成0的情况,虽然闭环系统一成之后这里的误差很难为0,但是实属是一个隐患.所以在这里我选择了在更新完x[n]e之后再移位.

再者,实际上,在后续的更新滤波器系数依然涉及到对误差的截位,也需要我们去操作一下.

编写testbeach..?

这里本来是不想写的…但是之前的一些博客被别人问了一下,这里就再写写怎么读入文件数据进行仿真.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
parameter 	data_num =500_000;
integer i;   //数组坐标
reg signed [11:0] stimulus[1:data_num];  //数组形式存储读出的数据
reg signed [11:0] misc_simulation[1:data_num];  //数组形式存储读出的数据
initial 
begin
   RST_N = 1'b1;  d_in = 12'b0; x_in = 12'b0;
	#60 RST_N = 1'b0;
	#60 RST_N = 1'b1; 
	$display("Running testbench");                       
    $readmemb("miscsig.txt", misc_simulation);  //将txt文件中的数据存储在数组中
	 $readmemb("xn.txt", stimulus);  //将txt文件中的数据存储在数组中

    i = 0;
    repeat(data_num) begin   //重复读取数组中的数据
        i = i + 1;
	    x_in= stimulus[i%data_num];
	    d_in = misc_simulation[i%data_num];
        #PERIOD;         //每个时钟读取一次
    end
	 $stop;
end

readmemb的意思是以二进制的形式读取一个数据,同样有readmemh等一系列函数.这里过于简单,不再赘述.

Matlab生成仿真数据

这里为了作死,我选取了一个双音信号做噪声输入,正常信号为一正弦波,还加入了两个相位差:
以下为系统参数,具体代码在之前的IIR滤波器的FPGA实现给出过.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 %=============设置系统参数==============%
f1=2000e3;        %设置波形频率
f2=100e3;
f3=10000e3
Fs=50e6;        %设置采样频率
L=50e3;         %数据长度
N=12;           %数据位宽
phase = pi/6;
%=============产生输入信号==============%
t=0:1/Fs:(1/Fs)*(L-1);
y1=sin(2*pi*f1*t);							%x[n]1
y2=1*sin(2*pi*f2*t+phase);     
y3=1*sin(2*pi*f3*t+phase);		%x[n]2
y4= y1+y2+y3;
y_n=round(y4*(2^(N-3)-1));      %N比特量化;如果有n个信号相加，则设置（N-n）

我们现在要从这个双音信号的干扰中还原出100Khz的这个信号出来.

Modelsim仿真结果

收敛时间

ubuntu截图
如果看系数的话,这大概在1ms左右的时间就收敛了,如果看波形的话,实际上大概在20us左右就收敛了.

滤波效果

说实话,滤波效果还不能说很好…
滤波效果
效果有点差强人意吧,不过大家其实也不必对四阶FIR滤波器抱太大的希望…

改进方向

实际上,在真正的电赛中,对频带有很严格的要求,但是收敛时间的要求也有点苛刻,但是也并不是真就做不了

此处只给方向,

阶数提高

这个针对的是频带,因为实际的滤波效果,滤波器性能和阶数是严格正相关的,所以提高阶数能把性能指标给提上去.但是可能会发生一方溢出,多方震荡的麻烦情节,增加许多Debug时间.

流水线化

这个针对的是系统的最高运行频率,但是流水线所带来的延时降低这个闭环系统的稳定性,所以这个流水性的阶数和步进因子u都需要好好地选择

变步长算法

这个针对的是收敛时间,通过判断误差的大小来选择不同的步长,从而实现更快的收敛速率.不过这个选择同样也需要小心翼翼地选择系数和阈值,否则可能会造成环路失控.

多相分解

这个针对的是收敛时间和滤波器性能.这个做法就是典型的”用空间换时间”,通过对输入信号的多相分解,实现高阶数的滤波效果,但是空间复杂度确实猛猛地涨.

结语

闭环系统的调试实在是太鸡儿难了,不过一开始还是因为自己理论还没过关,主要是对理论的指标还没过关,然后没有拆环分析,就一直瞎调,真就一顿操作猛如虎,一看波形丑如狗…

大二的时候都不知道自己怎么做出来的.突然想起了浣沙老师说过的: “你看上去他的代码写得是特别的烂,但是能用,以后工作了,也别想着重构这些代码,能用就行”.大二的那个版本做了变步长和多相分解,但是代码着实写得太烂了…并且用这个代码改上去也不难.我记得当时是根据两三篇硕士论文写出来的,有空再给大家贴上.

如果你觉得有丶收获的话