最小均方算法（使用 C++）。

• 下载源码 - 930 B

引言

由于其简单性，LMS 算法是最受欢迎的自适应算法之一。计算 LMS 不需要计算相关矩阵，甚至不需要计算矩阵逆。事实上，正是 LMS 算法的简单性使其成为衡量其他自适应滤波算法的标准。

最小均方算法的结构和操作概述

最小均方 (LMS) 算法是一种线性自适应滤波算法，由两个基本过程组成

1. 一个滤波过程，包括 (a) 计算由一组抽头输入产生的横向滤波器的输出，以及 (b) 通过将该输出与期望响应进行比较来生成估计误差。
2. 一个自适应过程，涉及根据估计误差自动调整滤波器的抽头权重。

自适应结构已用于自适应滤波中的不同应用，例如

• 噪声消除
• 系统辨识
• 自适应预测器
• 均衡
• 逆建模
• 干扰消除

在本文中，我们将实现我们在系统辨识中的示例。

系统辨识

图 1 显示了一种可用于系统辨识或建模的自适应滤波器结构。相同的输入 u 与自适应滤波器并行地输入到一个未知系统。误差信号 e 是未知系统 d 的响应与自适应滤波器 y 的响应之间的差值。误差信号反馈到自适应滤波器，并用于更新自适应滤波器的系数，直到总输出 y = d。当这种情况发生时，自适应过程结束，e 接近于零。

公式推导

正如您在图 1 中所看到的，y (n) 是自适应滤波器的输出

其中 wk (n) 代表特定时间 n 的 M 个权重或系数。卷积方程可以这样实现

for (int i = 0; i < M; i++)
        Y += (W[i] * X[i]);		//calculate filter output

需要一个性能指标来确定滤波器的优劣。该指标基于误差信号...

... 它是期望信号 d (n) 和自适应滤波器输出 y (n) 之间的差值。权重或系数 wk (n) 被调整以最小化均方误差函数。这个均方误差函数是 E [e2 (n)]，其中 E 代表期望值。由于有 k 个权重或系数，因此需要均方误差函数的梯度。可以使用 e2 (n) 的梯度来代替找到一个估计，得到

这可以这样实现

E = D - Y;		//calculate error signal

for (int i = 0; i < M ; i++)
	W[i] = W[i] + (mu * E * X[i]);		//update filter coefficients

代码中的 LMS 示例

我们使用图 1 中的自适应结构来说明自适应过程的以下步骤

1. 为 LMS 滤波器输入生成一些随机数据。
2. 假设一个我们要估计的系统，如下所示：H[M] = { 1, 0.5, 0.25, 0.125, 0.0625 }
3. 通过将生成的随机数据和假设的 H 进行卷积来构建期望信号。
4. 计算自适应滤波器输出 y。
5. 计算误差信号
6. 更新每个 LMS 滤波器权重
7. 为下一个输出采样点重复整个自适应过程。

生成输入和期望信号

for (int i = 0; i < I; i++)
		Input[i] = rand() / (float)RAND_MAX;

	for (int i = 0; i < I; i++)
		for (int j = 0; j < M; j++)
			if (i - j >= 0)
				Desired[i] += Input[i - j] * H[j];

总而言之

#define mu 0.2f				//convergence rate
#define M 5					//order of filter
#define I 1000				//number of samples

double Input[I] = { 0.0 };
double Desired[I] = { 0.0 };

//double H[M] = { 1, 0.5, 0.25, 0.125, 0.0625 };	//the main system
double H[M] = { 0.0625, 0.125, 0.25, 0.5, 1 };		//we need inverse of main system to convolution


void initialize()
{
	for (int i = 0; i < I; i++)
		Input[i] = rand() / (float)RAND_MAX;

	for (int i = 0; i < I; i++)
	for (int j = 0; j < M; j++)
	if (i - j >= 0)
		Desired[i] += Input[i - j] * H[j];
}

void main()
{
	initialize();

	long T, n = 0;
	double D, Y, E;
	double W[M] = { 0.0 };
	double X[M] = { 0.0 };

	FILE	*Y_out, *error, *weights;

	Y_out = fopen("Y_OUT", "w++");				//file for output samples
	error = fopen("ERROR", "w++");				//file for error samples
	weights = fopen("WEIGHTS", "w++");			//file for weights samples


	for (T = 0; T < I; T++)
	{
			for (int m = T; m > T - M; m--){
				if (m >= 0)
					X[M + (m - T) - 1] = Input[m];	//X new input sample for 
									//LMS filter
				else break;
			}

			D = Desired[T];					//desired signal
			Y = 0;						//filter’output set to zero

			for (int i = 0; i < M; i++)
				Y += (W[i] * X[i]);			//calculate filter output

			E = D - Y;					//calculate error signal

			for (int i = 0; i < M; i++)
				W[i] = W[i] + (mu * E * X[i]);		//update filter coefficients

			fprintf(Y_out, "\n % 10g % 10f", (float)T, Y);
			fprintf(error, "\n % 10g % 10f", (float)T, E);
	}

	for (int i = 0; i < M; i++)
		fprintf(weights, "\n % 10d % 10f", i, W[i]);

	fclose(Y_out);
	fclose(error);
	fclose(weights);
}

参考文献

• 自适应滤波器，Simon Haykin
• 使用 C 和 TMS320C6x DSK 的 DSP 应用