Administrator
发布于 2024-06-08 / 0 阅读
0
0

应用随机梯度算法解决信道均衡问题

简介

利用优化算法并构建滤波器来解决信号通道的失真(噪声、多径传播)等问题。在数字通信中,发射机向接收机传输一个符号序列 {s(n)}_{n>0},这些符号取值为 1-1,来编码有用的信息。但是在传输过程中,通常信号会失真,这种失真通常可被视为,与传播信道对应滤波器的卷积。

因此,在信号接收端,需要对该传播信道进行去卷积以恢复被传输的符号。但是在去卷积之前,必须对信道进行估计,也就是了解信道的具体特性(例如信道的冲激响应或频率响应)。为了了解特性,我们需要使用假设已知的训练序列。

训练序列是 发射机和接收机之间 预先约定的已知符号序列。接收端在收到经过信道传输的训练序列后,可以将其与原始序列进行比较,从而推断信道的特性。这些信道特性我们称之为“信道系数” w(k)。在信道系数估计完成后,我们就可以设计滤波器来实现去卷积操作,来减小由于信道引起的失真。

在这种情况下,随机梯度算法(LMS)被经典地用于该方面。因此,后续我们将使用 LMS 应用于一个简化的信道均衡问题。

传输链路的建模与仿真

w(k) = \begin{cases} \frac{1}{2} \left( 1 + \cos\left(\frac{2\pi}{\beta}(k - 1)\right) \right), & \text{for } k = 0, 1, 2, \\ 0, & \text{else} \end{cases}

此外,接收信号受到加性高斯噪声 u(n) 的干扰,加性高斯噪声具有白噪音的特性(独立同分布),零均值,方差为 \sigma^2 = 0.001

因此,信号的表达式为:

x(n) = \sum_{k=0}^{2} w(k)s(n-k) + u(n)

其中, \beta 控制信道引入的失真水平。

我们先在 Matlab 中使用函数 stem 生成由符号 1-1 组成的随机序列信号 s(n),并假设 1-1 取值是等概率的。我们可以使用函数 randn,并假设信号的长度为 N = 1500

clear all
clc
close all

% 参数设置
N = 1500;

% 生成随机信号 s(n) 为 +1 或 -1
s = sign(randn(N,1));

% 绘制离散信号 s(n)
figure;
stem(s); % 离散信号绘图
title('信号 s(n)', 'FontSize', 16);
xlabel('信号长度 N', 'FontSize', 14);
ylabel('s(n)', 'FontSize', 14);
set(gca, 'FontSize', 12); % 放大坐标轴刻度

Alt Text

只展示了 200 的序列

由于信号 s(n) 是由符号 1-1 组成的随机序列,且被假设等概率的:

X = {+1, -1}

p_{+1} = p_{-1} = \frac{1}{2}

因此,s(n) 的期望形式为:

E[s(n)] = p_{+1} \cdot (+1) + p_{-1} \cdot (-1)
E[s(n)] = \frac{1}{2} \cdot (+1) + \frac{1}{2} \cdot (-1)
E[s(n)] = 0

我们接下来绘制传播信道,其公式为:

w(k) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{\beta}(k - 1) \right) \right)

\beta 取不同值: 0.25, 24 来进行观察

对于 \beta = 0.25

\beta 较小时,例如 0.25,传输信道具有较短的脉冲响应,这会导致余弦函数的快速变化。在这种情况下,可以预期 w(k)k = 0k = 2 之间有显著变化。

w(0) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{0.25}(-1) \right) \right) = \frac{1}{2} \left( 1 + \cos(-8\pi) \right) = \frac{1}{2}(1 + 1) = 1
w(1) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{0.25}(0) \right) \right) = \frac{1}{2} \left( 1 + \cos(0) \right) = \frac{1}{2}(1 + 1) = 1
w(2) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{0.25}(1) \right) \right) = \frac{1}{2} \left( 1 + \cos(8\pi) \right) = \frac{1}{2}(1 + 1) = 1

因此,当 \beta = 0.25 时,w(k) 在点 k = 0, 1, 2 处的值均为 1

Alt Text

可以看到传播信道在前 3 个采样点上具有恒定的脉冲响应。

对于 \beta = 2

w(0) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{2}(-1) \right) \right) = \frac{1}{2} \left( 1 + \cos(-\pi) \right) = \frac{1}{2}(1 - 1) = 0
w(1) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{2}(0) \right) \right) = \frac{1}{2} \left( 1 + \cos(0) \right) = \frac{1}{2}(1 + 1) = 1
w(2) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{2}(1) \right) \right) = \frac{1}{2} \left( 1 + \cos(\pi) \right) = \frac{1}{2}(1 - 1) = 0

因此,对于 \beta = 2w(k) 在点 k = 0k = 2 的值为 0,而在点 k = 1 的值为 1

Alt Text

可以看到传播信道的脉冲响应主要集中在 k = 1,主要影响第一个延迟的信号。

对于 \beta = 4

\beta 更大时,例如 4,传输信道的脉冲响应周期进一步延长,余弦函数的变化变得更慢。在这种情况下,可以预期 w(k)k = 0k = 2 之间的变化较小。

w(0) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{4}(-1) \right) \right) = \frac{1}{2} \left( 1 + \cos\left(-\frac{\pi}{2}\right) \right) = \frac{1}{2}(1 + 0) = \frac{1}{2}
w(1) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{4}(0) \right) \right) = \frac{1}{2} \left( 1 + \cos(0) \right) = \frac{1}{2}(1 + 1) = 1
w(2) = \frac{1}{2} \left( 1 + \cos \left( \frac{2\pi}{4}(1) \right) \right) = \frac{1}{2} \left( 1 + \cos\left(\frac{\pi}{2}\right) \right) = \frac{1}{2}(1 + 0) = \frac{1}{2}

因此,对于 \beta = 4w(k) 在点 k = 0k = 2 的值为 0.5,而在点 k = 1 的值为 1

Alt Text

可以看到传播信道具有一个对称的脉冲响应,其峰值位于 k = 1,如 \beta = 2 的情况(图 3),但它在第一个和第三个采样点上也包含非零值。即脉冲响应更加分散且对称,但主要集中在中心。

这表明,具有均匀或扩展脉冲响应的信道需要均衡器来补偿分布在多个采样点上的失真;而具有集中脉冲响应的信道可能更容易均衡,因为失真主要局限于单个采样点。

后续中,我们将使用 \beta = 0.25,通过已经构建好的随机序列信号 s(n),使用 Matlab 函数 filter 来添加加性噪声,以获得信号 x(n)

%% 给信号加噪音
% 滤波器脉冲响应(权值)
w = [1 1 1];

% 对 s(n) 进行滤波
% filter(num, den, data_to_filter)
x = filter(w, 1, s);

sigma_u = sqrt(0.001); % 高斯噪声标准差

% 向滤波后的信号中加入高斯噪声
u = sigma_u * randn(N,1);
x = x + u;

% 绘制加入噪声后的信号 x(n)
figure;
stem(x);
title('信号 x(n)', 'FontSize', 16);
xlabel('信号长度 N', 'FontSize', 14);
ylabel('x(n)', 'FontSize', 14);
set(gca, 'FontSize', 12); % 放大坐标轴刻度

Alt Text

实现均衡

为了实现信道均衡,我们最终目标是确定最优维纳滤波器 h_{\text{opt}} 的理论表达式,以最小化均方误差。为此,首先需要计算信号 x(n) 的自相关函数和 x(n)d(n) 的互相关函数,分别记为 r_{xx}(k)r_{dx}(k)。对于真实且平稳的信号,其相关函数公式为:

r_{xx}(k) = E[x(n)x(n-k)]
r_{dx}(k) = E[d(n)x(n-k)]

我们首先研究(由传播信道卷积得到的)发射信号 m(n) 的自相关函数 r_{mm}(k),其中 m(n) = \sum_{k=0}^{2} w(k)s(n-k)。首先计算 r_{mm}(0)r_{mm}(1)r_{mm}(2) 来观察结果。同时注意,我们之前假设符号 s(n) 为相互独立的。

计算 r_{mm}(0)

r_{mm}(0) = E \left[ \left( \sum_{k=0}^{2} w(k)s(n-k) \right)^2 \right]
r_{mm}(0) = \sum_{k=0}^{2} w^2(k) E[s^2(n-k)]

由于 s(n) 是相互独立的,因此有:

E[s^2(n-k)] = 1

由于 s(n) 的特性:

E[s^2(n)] = 1 \cdot P(s(n) = 1) + 1 \cdot P(s(n) = -1)
E[s^2(n)] = 1 \cdot \frac{1}{2} + 1 \cdot \frac{1}{2} = 1

因此:

r_{mm}(0) = \sum_{k=0}^{2} w^2(k)
r_{mm}(0) = w(0)^2 + w(1)^2 + w(2)^2 = 3

计算 r_{mm}(1)

r_{mm}(1) = E[m(n)m(n-1)]
r_{mm}(1) = E \left[ \left( \sum_{k=0}^{2} w(k)s(n-k) \right) \left( \sum_{j=0}^{2} w(j)s(n-1-j) \right) \right]
r_{mm}(1) = \sum_{k=0}^{2} \sum_{j=0}^{2} w(k)w(j) E[s(n-k)s(n-1-j)]

同样,仅当 k = j + 1 时,相关函数才不为零,因此:

r_{mm}(1) = w(0)w(1)E[s(n-1)s(n-1)] + w(1)w(2)E[s(n-2)s(n-2)]

由于 E[s(n)] = 0,因此:

E[s^2(n)] = 1

因此:

r_{mm}(1) = w(0)w(1) \cdot 1 + w(1)w(2) \cdot 1
r_{mm}(1) = w(0)w(1) + w(1)w(2) = 2

计算 r_{mm}(2)

r_{mm}(2) = E \left[ \left( \sum_{k=0}^{2} w(k)s(n-k) \right) \left( \sum_{j=0}^{2} w(j)s(n-2-j) \right) \right]
r_{mm}(2) = \sum_{k=0}^{2} \sum_{j=0}^{2} w(k)w(j) E[s(n-k)s(n-2-j)]

同样,仅当 k = j + 2 时,相关函数才不为零,因此:

r_{mm}(2) = w(0)w(2)E[s(n-2)s(n-2)]

由于 E[s(n)] = 0,因此:

E[s^2(n-2)] = 1

因此:

r_{mm}(2) = w(0)w(2) = 1

接下来我们需要找到当 k > 2 时,r_{mm}(k) = 0 并推出任意 kr_{mm}(k) 的表达式。

为了验证 r_{mm}(k) = 0k > 2,我们先计算 r_{mm}(3) 看看:

r_{mm}(3) = E \left[ \left( \sum_{k=0}^{2} w(k)s(n-k) \right) \left( \sum_{j=0}^{2} w(j)s(n-3-j) \right) \right]
r_{mm}(3) = \sum_{k=0}^{2} \sum_{j=0}^{2} w(k)w(j) E[s(n-k)s(n-3-j)]

同样,仅当 k = j + 3 时,相关函数才不为零,因此:

r_{mm}(3) = w(0)w(3)E[s(n-3)s(n-3)]

由于:

w(k) = 0 , \text{对于 } k \neq 0, 1, 2

因此:

w(3) = 0

因此:

r_{mm}(3) = 0

结论

r_{mm}(k) = \begin{cases} w(0)^2 + w(1)^2 + w(2)^2 & \text{当 } k = 0 \\ w(0)w(1) + w(1)w(2) & \text{当 } k = 1 \\ w(0)w(2) & \text{当 } k = 2 \\ 0 & \text{当 } k > 2 \end{cases}
r_{mm}(k) = \begin{cases} 3 & \text{当 } k = 0 \\ 2 & \text{当 } k = 1 \\ 1 & \text{当 } k = 2 \\ 0 & \text{当 } k > 2 \end{cases}

计算 r_{xx}(k)

我们使用已经得到的结果来计算 r_{xx}(k)

我们已知 x(n) = m(n) + u(n),并且噪声 u(n) 被假设为独立于传输符号的。

r_{xx}(k) = E[x(n)x(n-k)]
r_{xx}(k) = E\left[(m(n) + u(n))(m(n-k) + u(n-k))\right]
r_{xx}(k) = E[m(n)m(n-k)] + E[m(n)u(n-k)] + E[u(n)m(n-k)] + E[u(n)u(n-k)]

由于 m(n)u(n) 是独立的:

E[m(n)u(n-k)] = 0

并且 u(n) 是零均值的白噪声:

E[u(n)u(n-k)] = \sigma_u^2\delta(k)

因此:

r_{xx}(k) = r_{mm}(k) + \sigma_u^2\delta(k)
r_{xx}(k) = \begin{cases} w(0)^2 + w(1)^2 + w(2)^2 + \sigma_u^2 & \text{当 } k = 0 \\ w(0)w(1) + w(1)w(2) & \text{当 } k = 1 \\ w(0)w(2) & \text{当 } k = 2 \\ 0 & \text{当 } k > 2 \end{cases}
r_{xx}(k) = \begin{cases} 3 + \sigma_u^2 & \text{当 } k = 0 \\ 2 & \text{当 } k = 1 \\ 1 & \text{当 } k = 2 \\ 0 & \text{当 } k > 2 \end{cases}

h_{\text{opt}} 的表达式

最终我们计算 r_{dx}(k) 并最终给出 h_{\text{opt}} 的表达式。

r_{dx}(k) = E[d(n)x(n-k)]

已知 d(n) = s(n-2)x(n) = m(n) + u(n)

r_{dx}(k) = E[s(n-2)(m(n-k) + u(n-k))]

由于 u(n)s(n) 是独立的且 u(n) 的均值为 0

E[s(n-2)u(n-k)] = 0

因此:

r_{dx}(k) = E[s(n-2)m(n-k)]

由于已知:

m(n) = \sum_{k=0}^{2} w(k)s(n-k)

因此:

r_{dx}(k) = E[s(n-2)(w(0)s(n-k) + w(1)s(n-k-1) + w(2)s(n-k-2))]

k = 0 时:

r_{dx}(0) = w(2)E[s(n-2)s(n-2)] = w(2)

k = 1 时:

r_{dx}(1) = w(1)E[s(n-2)s(n-2)] = w(1)

k = 2 时:

r_{dx}(2) = w(0)E[s(n-2)s(n-2)] = w(0)

结论

r_{dx}(k) = \begin{cases} w(2) & \text{当 } k = 0 \\ w(1) & \text{当 } k = 1 \\ w(0) & \text{当 } k = 2 \\ 0 & \text{当 } k > 2 \end{cases}

并且:

h_{\text{opt}} = R_{xx}^{-1}r_{dx}
h_{\text{opt}} = \begin{pmatrix} r_{xx}(0) & r_{xx}(1) & r_{xx}(2) \\ r_{xx}(1) & r_{xx}(0) & r_{xx}(1) \\ r_{xx}(2) & r_{xx}(1) & r_{xx}(0) \end{pmatrix}^{-1} \begin{pmatrix} r_{dx}(0) \\ r_{dx}(1) \\ r_{dx}(2) \end{pmatrix}
h_{\text{opt}} = \begin{pmatrix} r_{xx}(0) & r_{xx}(1) & r_{xx}(2) \\ r_{xx}(1) & r_{xx}(0) & r_{xx}(1) \\ r_{xx}(2) & r_{xx}(1) & r_{xx}(0) \end{pmatrix}^{-1} \begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}

实现

我们将在 Matlab 中设计一个滤波器,并使用两种方法计算自相关和互相关来构建滤波器:一种方法是基于理论公式,用 Toeplitz 函数构造自相关矩阵;另一种方法是通过 Matlab 的 xcorr 函数以数值方式计算这些相关函数,然后比较两种方法下得到的滤波器结果。

我们先看第一种实现,使用上述得到的计算公式,并利用 Toeplitz 函数构造自相关矩阵。

%% 第一种实现
L = 11;

% 计算自相关函数 r_xx(k)
% r_xx(0)
r_xx1 = w(1)^2 + w(2)^2 + w(3)^2 + sigma_u^2;
% r_xx(1)
r_xx2 = w(1)*w(2) + w(2)*w(3);
% r_xx(2)
r_xx3 = w(1)*w(3);
% 对 k > 2,r_xx(k) = 0
r_xx = [r_xx1; r_xx2; r_xx3; zeros(8,1)];

% 构建 Toeplitz 矩阵 R_xx
R_xx = toeplitz(r_xx);

% 计算互相关函数 r_dx(k)
% r_dx(0)
r_dx1 = w(3);
% r_dx(1)
r_dx2 = w(2);
% r_dx(2)
r_dx3 = w(1);
% 对 k > 2,r_dx(k) = 0
r_dx = [r_dx1; r_dx2; r_dx3; zeros(8,1)];

% 使用 Wiener-Hopf 方程求解最优滤波器系数
h_opt = inv(R_xx)*r_dx;
disp('Wiener 最优滤波器 h_opt:');
disp(h_opt);


第二种实现,通过 Matlab 的 xcorr 函数以数值方式(不依赖理论公式,直接应用样本数据)计算这些相关函数

%% 第二种实现
hopt_theo = h_opt;

% 定义 d(n) = s(n-2)
d = s(1:end-2);
x2 = x(3:end);

% 数值计算 r_xx
r_xx_num = xcorr(x2);
r_xx_num = r_xx_num(N-2:N-2+L-1);

% 数值计算 r_dx
r_dx_num = xcorr(d, x2);
r_dx_num = r_dx_num(N-2:N-2+L-1);

% 构建数值 Toeplitz 矩阵 R_xx_num
R_xx_num = toeplitz(r_xx_num);

% 数值求解最优滤波器系数
h_opt_num = inv(R_xx_num)*r_dx_num;

% 绘制理论与数值求解得到的最优滤波器系数
figure,
stem(hopt_theo), title('理论最优滤波器', 'FontSize',16), hold on;
stem(h_opt_num), title('数值最优滤波器', 'FontSize', 16)
xlabel('滤波器阶数 L','FontSize',14),
ylabel('幅度', 'FontSize',14)
set(gca, 'FontSize', 14);
legend('理论结果','数值结果')

Alt Text

我们发现数值滤波器(Numérique)非常接近理论最优滤波器(Théorique)。

现在,我们应用 LMS 算法来迭代解决信道估计问题,这里假设 L = 11。我们创建一个函数 algoLMS,该函数接收输入信号 x(n)d(n),脉冲响应长度 L,以及步长 \mu,并输出后验误差序列 e^+(n) 和通过 LMS 算法迭代得到的脉冲响应 h_n。并通过调整步长 \mu 来测试算法,并且展示滤波器系数误差的范数 |h_n - h_{\text{opt}}|_2

%% LMS算法
% 学习步长数组
mu = [0.001, 0.005, 0.01, 0.02];

% 对不同步长 mu 执行 LMS 算法
[err1, h1] = algoLMS(x2, d, L, mu(1));
[err2, h2] = algoLMS(x2, d, L, mu(2));
[err3, h3] = algoLMS(x2, d, L, mu(3));
[err4, h4] = algoLMS(x2, d, L, mu(4));

% 计算滤波器系数误差的范数
normeh1 = sqrt(sum((h1 - hopt_theo).^2));
normeh2 = sqrt(sum((h2 - hopt_theo).^2));
normeh3 = sqrt(sum((h3 - hopt_theo).^2));
normeh4 = sqrt(sum((h4 - h_opt).^2, 1));

% 绘制后验误差 e+(n) 随 n 的变化
% 使用子图分别展示不同 μ 值下的后验误差,以获得更清晰的对比

figure;

% μ = 0.001
subplot(2,2,1);
plot(err1); grid on;
xlabel('信号长度 N', 'FontSize', 14);
ylabel('e(n)', 'FontSize', 14);
title('μ = 0.001', 'FontSize', 14);

% μ = 0.005
subplot(2,2,2);
plot(err2); grid on;
xlabel('信号长度 N', 'FontSize', 14);
ylabel('e(n)', 'FontSize', 14);
title('μ = 0.005', 'FontSize', 14);

% μ = 0.01
subplot(2,2,3);
plot(err3); grid on;
xlabel('信号长度 N', 'FontSize', 14);
ylabel('e(n)', 'FontSize', 14);
title('μ = 0.01', 'FontSize', 14);

% μ = 0.02
subplot(2,2,4);
plot(err4); grid on;
xlabel('信号长度 N', 'FontSize', 14);
ylabel('e(n)', 'FontSize', 14);
title('μ = 0.02', 'FontSize', 14);


% 绘制与 h_opt 的误差范数随 n 的变化
figure;
plot(normeh1); grid on; hold on;
plot(normeh2);
plot(normeh3);
plot(normeh4);
xlabel('信号长度 N', 'FontSize', 14);
ylabel('||h(n)-h_{opt}||^2', 'FontSize', 14);
title('后验误差范数', 'Interpreter', 'none', 'FontSize', 16);
legend('\mu = 0.001', '\mu = 0.005', '\mu = 0.01', '\mu = 0.02');

子函数:

function [eplus,h_tout] = algoLMS(x,d,L,mu)
% 实现 LMS 算法的函数

N = length(x);
h = ones(L,1);
eplus = zeros(N,1);
h_tout = zeros(L,N);

for n = L:N
    xn = x(n:-1:n-L+1);
    yn = h'*xn;
    eplus(n) = d(n) - yn;
    % LMS 更新方程
    h = h + mu*xn*(d(n)-xn'*h);
    h_tout(:,n) = h;
end

end

Alt Text

Alt Text

我们下面计算理论上使算法收敛的步长 \mu 的最大值,并和实际值做对比。

收敛的充分必要条件为:

0 < \mu < \frac{2}{\lambda_{\text{max}}}

其中, \lambda_{\text{max}} 是自相关矩阵 R_{xx} 的最大特征值。

我们有:

E[h_n] = E[h_{n-1}] - \mu E[x(n)x(n)^T h_{n-1}] + \mu E[d(n)x(n)]

在假设 x(n)h_{n-1} 独立的情况下,可以写为:

E[x(n)x(n)^T h_{n-1}] = E[x(n)x(n)^T]E[h_{n-1}]

因此:

E[h_n] = E[h_{n-1}] - \mu E[x(n)x(n)^T]E[h_{n-1}] + \mu E[d(n)x(n)]

此外:

\mu E[d(n)x(n)] = R_{xx} \cdot h_{\text{opt}}

因此:

E[h_n] - h_{\text{opt}} = (I - \mu R_{xx})(E[h_{n-1}] - h_{\text{opt}})

假设:

\delta h_n = E[h_n] - h_{\text{opt}}

因此:

\delta h_n = (I - \mu R_{xx})\delta h_{n-1}

经过对 R_{xx} 的一系列变换,我们得到:

\delta h_n = (I - \mu R_{xx})\delta h_{n-1}
\delta h_n(l) = (1 - \mu \lambda_l)^n \delta h_0(l)
|\mu \lambda_l| < 1 \iff 0 < \mu < \frac{2}{\lambda_l}

结论:

0 < \mu < \frac{2}{\lambda_{\text{max}}}

理论最大值为:

\mu = \frac{2}{\lambda_{\text{max}}} = \frac{2}{8.6236} = 0.2319

其中 \lambda_{\text{max}} 由 Matlab 确定。详细解释后续补充

% 计算 R_xx 的特征值
lambda = eig(R_xx);

% 找到最大特征值
lambda_max = max(lambda);

% 根据公式计算理论上使算法收敛的最大步长 mu_theorique_max
mu_theorique_max = 2 / lambda_max;

% 显示结果
disp('理论最大步长 mu_theorique_max:');
disp(mu_theorique_max);

由于理论计算是在以下假设条件下完成的:输入信号是平稳的,信道模型是精确的。而在实践中,真实信号并不总是满足这些假设。因此,应适当降低 \mu 的值。在我们的案例中,选择 \mu = 0.01

我们现在改变滤波器的阶数 L 和噪声水平 \sigma^2_u,以测试算法的鲁棒性。特别是,绘制误差范数 |h_n - h_{\text{opt}}|_2\sigma^2_u 的变化曲线。后续待补充。


评论