本帖最后由 wopt 于 2015-4-9 21:16 编辑
3.7 LMS
Least mean-square algorithm，最小均方算法。由Widrow等人在1975年提出。
定义：
u(n)：输入向量，或称为训练样本
w(n)：权值向量
e(n)：偏差
d(n)：期望响应
d^(n)：对期望响应d(n)的估计
n：迭代次数
µ：步长因子或步长参数。
LMS算法：
步骤1：
初始化
迭代次数：n = 0
权向量：w(0) = 0
估计误差：e(0) = d(0) - d^(0) = d(0)
输入向量：u(0) = [u(0) 0 …… 0]T
步骤2：对 n = 0, 1 , 2…
权向量的更新：w(n + 1) = w(n) + µ * u(n) * e*(n)
期望信号的估计：d^(n + 1) = w(n + 1) * u(n + 1)
估计误差：e(n + 1) = d(n + 1) - d^(n + 1)
步骤3： n = n + 1，重复步骤2
先来一个简单的例子，熟悉下lms算法。Matlab上菜。
clc; clear all;
%% lms
N = 2048; % 信号的取样点数
M = 10; % 抽头的个数
iter = N; % 迭代次数
t = 1 : N;
s = 5 * sin(2 * pi * 0.01 * t); % 源信号
x = s + 0.5 * randn(1, N); % 输入信号
x = x';
d = s'; % 期望信号
e = zeros(iter, 1);
w = zeros(M, iter);
rho_max = max(eig(x * x')); % 输入信号相关矩阵的最大特征值
mu = rand(1) * (2 / rho_max) ; % 收敛因子 0 < mu < 2/rho
for n = M + 1 : iter -1
w( :, n + 1) = w( :, n) + mu * x(n - 1 : -1 : n - M) * conj(e(n));
y(n + 1) = w( :, n + 1)' * x(n : -1 : n - M + 1);
e(n + 1) = d(n + 1, : ) - y(n + 1);
end
%% 画输入输出信号以及输出信号和预期输出信号的误差
figure(1);
subplot(2, 1, 1);
plot(t, x);
grid;
title('输入信号');
xlabel('样点数');
ylabel('幅值');
subplot(2,1,2);
plot(t, y);
grid;
title('输出信号');
xlabel('样点数');
ylabel('幅值');
% 绘制自适应滤波器输出信号,预期输出信号和两者的误差
figure(2)
subplot(3, 1, 1);
plot(t, d);
grid;
title('预期输出');
xlabel('样点数');
ylabel('幅值');
subplot(3, 1, 2);
plot(t, y);
grid;
title('自适应滤波器输出');
xlabel('样点数');
ylabel('幅值');
subplot(3, 1, 3);
plot(t, e);
grid;
title('误差');
xlabel('样点数');
ylabel('幅值');
在来一个难一点的。引用《现代数字信号处理及其应用》—何子述、夏威 课后的一个仿真题。
clc; clear all;
%% 产生512点输入序列
M = 2; % AR 阶数
N = 512; % 样本数
trials = 100; % 试验次数
a1 = -0.975; % AR parameter
a2 = 0.95;
var = 0.0731; % 噪声方差
v = 0 + sqrt(var).* randn(N, 1); % 产生v(n)
u0 = [0, 0]; % u(0)初始化
a = [1 a1 a2];
b = 1;
zi = filtic(b, a, u0); % 初始化滤波器 直接II型
[u, zf] = filter(b, a, v, zi); % 产生观测信号u(n)
%% 用LMS算法计算w1和w2
mu1 = 0.05;
mu2 = 0.005;
w1 = zeros(M, N); % 初始化权向量
w2 = zeros(M, N);
w1_trials = zeros(M, N, trials);
w2_trials = zeros(M, N, trials);
e1 = zeros(N, 1); % 初始化误差
e2 = zeros(N, 1);
e1_trials = zeros(N, 1, trials);
e2_trials = zeros(N, 1, trials);
d1 = zeros(N, 1); % 初始化期望响应
d2 = zeros(N, 1);
wopt = zeros(2, trials);
Jmin = zeros(1, trials);
sum_eig = zeros(1, trials);
for i = 1 : trials
v = 0 + sqrt(var).* randn(N, 1);
zi = filtic(b, a, u0);
[u, zf] = filter(b, a, v, zi);
uu = zeros(M, N - M + 1);
for k = 1 : (N - M + 1)
uu( :, k) = u(k + M - 1 : -1 : k).';
end
R = uu * uu' / (N - M + 1); % u(n)自相关矩阵
p = sum(uu(:, 1 : N - M) * diag(u(3 : N)), 2) / (N - M);
[v, d] = eig(R);
wopt( :, i) = inv(R) * p;
Jmin(i) = R(1, 1) - p' * wopt( :, i); % Jmin
sum_eig(i) = d(1, 1) + d(2, 2); % 特征值之和
for n = 3 : N -1
w1( :, n + 1) = w1( :, n) + mu1 * u(n - 1 : -1 : n - 2, :) * conj(e1(n));
d1(n + 1) = w1( :, n + 1)' * u(n : -1 : n - 1, : );
e1(n + 1) = u(n + 1, : ) - d1(n + 1);
w2( :, n + 1) = w2( :, n) + mu2 * u(n - 1 : -1 : n - 2, :) * conj(e2(n));
d2(n + 1) = w2( :, n + 1)' * u(n : -1 : n - 1, : );
e2(n + 1) = u(n + 1, : ) - d2(n + 1);
end
w1_trials ( :, :, i) = w1;
e1_trials ( :, :, i) = e1 .* conj(e1);
w2_trials ( :, :, i) = w2;
e2_trials ( :, :, i) = e2 .* conj(e2);
end
w1_trials_ave = sum(w1_trials , 3) ./ trials;
e1_trials_ave = sum(e1_trials , 3) ./ trials;
w2_trials_ave = sum(w2_trials , 3) ./ trials;
e2_trials_ave = sum(e2_trials , 3) ./ trials;
sJmin = sum(Jmin) / trials; % 维纳误差 平均
Jex1 = e1_trials_ave - sJmin; % 剩余均方误差
Jex2 = e2_trials_ave - sJmin; % 剩余均方误差
sum_eig_trials = sum(sum_eig) / trials; % 平均特征值之和
Jexfin1 = mu1 * sJmin * (sum_eig_trials / (2 - mu1 * sum_eig_trials));
Jexfin2 = mu2 * sJmin * (sum_eig_trials / (2 - mu2 * sum_eig_trials));
M1 = Jexfin1 / sJmin; % 失调参数
M2 = Jexfin2 / sJmin; % 失调参数
disp(M1);
disp(M2);
%% 画滤波器的权向量和学习曲线
figure(1); % 画权值迭代过程
plot(w1(1, :), 'b'); % 单次估计 w1
hold on
plot(w1(2, :), 'r');
hold on
plot(w1_trials_ave(1, :), 'b'); % 估计平均 w1
hold on
plot(w1_trials_ave(2, :), 'r');
xlabel('迭代次数n');
ylabel('抽头权值');
title('滤波器权系数的更新过程 \mu=0.05');
figure(2); % 画权值迭代过程
plot(w2(1, :), 'b'); % 单次估计 w2
hold on
plot(w2(2, :), 'r');
hold on
plot(w2_trials_ave(1, :), 'b'); % 估计平均 w2
hold on
plot(w2_trials_ave(2, :), 'r');
xlabel('迭代次数n');
ylabel('抽头权值');
title('滤波器权系数的更新过程 \mu=0.005');
figure(3); % 画学习曲线
plot(Jex1, 'b'); % 画mu1时的学习曲线
hold on
plot(Jex2, 'r'); % 画mu2时的学习曲线
hold on
axis([0,550,0,1]);
xlabel('迭代次数n');
ylabel('MSE');
title('LMS算法的学习曲线');
legend('\mu1=0.05', '\mu2=0.005');
接下来在看看开发板附赠程序LMS的代码#include"math.h"
#define PI 3.1415926
#define COEFFNUMBER 16
#define INPUTNUMBER 1024
int FIRLMS(int *nx,float *nh,int nError,int nCoeffNumber);
float h[COEFFNUMBER],fU;
int xx[INPUTNUMBER],rr[INPUTNUMBER],wc[INPUTNUMBER];
main()
{
int i,nLastOutput;
nLastOutput=0;
fU=0.0005;
for ( i=0;i<COEFFNUMBER;i++ ) h[i]=0;
for ( i=0;i<INPUTNUMBER;i++ )
{
xx[i]=256*sin(i*2*PI/34);
rr[i]=wc[i]=0;
}
for ( i=COEFFNUMBER+1;i<INPUTNUMBER;i++ )
{
nLastOutput=FIRLMS(xx+i,h,nLastOutput-xx[i-1],COEFFNUMBER); // break point
rr[i]=nLastOutput;
wc[i]=rr[i]-xx[i];
}
exit(0);
}
int FIRLMS(int *nx,float *nh,int nError,int nCoeffNumber)
{
int i,r;
float fWork;
r=0;
for ( i=0;i<nCoeffNumber;i++ )
{
fWork=nx[i]*nError*fU;
nh[i]+=fWork;
r+=(nx[i-i]*nh[i]);
}
r/=128;
return r;
}
// LMS FIR (delayed version)
ushort dlms (DATA *x, DATA *h, DATA *r, DATA *des,DATA *dbuffer, DATA step, ushort nh, ushort nx)
// Adaptive delayed LMS filter (fast implemented)
ushort oflag = dlmsfast (DATA *x, DATA *h, DATA *r, DATA*des, DATA *dbuffer, DATA step, ushort nh, ushort nx)
|
楼主
