1️⃣ LS（最小二乘）信道估计

OFDM系统的信道估计常在频域进行，因为OFDM本身就是基于频域的。频域模型可以表示为：
$$Y(f)=X(f)H(f)+Z(f)$$
其中，$Y(f)$表示接收信号的频域表示，$X(f)$表示发送信号的频域表示，$H(f)$表示系统的频率响应【对于OFDM系统，系统的频率响应是是单位脉冲响应的h(n)的DTFT】。为简化描述，下面的$Y(f)$，$X(f)$和$H(f)$分别表示为$Y$，$X$，$H$。

LS信道估计的目标是最小化下面的代价函数：
$$J({\hat{H}}_{LS})=\Vert Y-X{\hat{H}}_{LS}{\Vert }^2$$
其中，${\hat{H}}_{LS}$表示LS信道估计值

令上面的代价函数关于 $\hat{H}$ 的偏导数等于 0 ，即

$$\frac{\partial J({\hat{H}}_{LS})}{\partial {\hat{H}}_{LS}}=-2{\left(X^{\mathrm{H}}Y\right)}^{\ast }+2{\left(X^{\mathrm{H}}X\hat{H}\right)}^{\ast }=0$$

然后可以得到 $X^{H}X\hat{H}=X^{\mathrm{H}}Y$ ，由此得到LS信道估计的解为：

$${\hat{H}}_{\mathrm{LS}}={\left(X^{\mathrm{H}}X\right)}^{-1}{X}^{\mathrm{H}}Y=X^{-1}Y$$

在OFDM系统中，可以用导频信号来得到${\hat{H}}_{\mathrm{LS}}$。对于上式，X是发送的导频信号，Y是接收的导频信号。

由于导频符号的数量有限，通常需要对LS估计结果进行插值来估算完整的数据符号【要传输的OFDM符号】的信道估计

LS 信道估计的均方误差（MSE）为

$$\begin{array}{rl}{\mathrm{MSE}}_{\mathrm{LS}}& =E\left\{{\left(H-{\hat{H}}_{\mathrm{LS}}\right)}^{\mathrm{H}}\left(H-H_{\mathrm{LS}}\right)\right\}\\ & =E\left\{{\left(H-X^{-1}Y\right)}^{\mathrm{H}}\left(H-X^{-1}Y\right)\right\}\\ & =E\left\{{\left(X^{-1}Z\right)}^{\mathrm{H}}\left(X^{-1}Z\right)\right\}\\ & =E\left\{Z^{\mathrm{H}}{\left(X{X}^{\mathrm{H}}\right)}^{-1}Z\right\}\\ & =\frac{{\sigma }_z^2}{{\sigma }_x^2}\end{array}$$

注意，上式中的 MSE 与信噪比 ${\sigma }_x^2/{\sigma }_z^2$ 成反比，噪声越大，误差越小。这意味着LS估计增强了噪声，在信道处于深度衰落时更是如此。然而，LS 方法由于简单而被广泛应用于信道估计。

---

2️⃣ MMSE（最小均方误差）信道估计

MMSE估计是LS估计的改进，在LS估计的基础上乘以加权矩阵$W$，使之与实际信道参数更为接近。

同样的，对于OFDM系统，频域模型可以表示为：
$$Y(f)=X(f)H(f)+Z(f)$$
其中，$Y(f)$表示接收信号的频域表示，$X(f)$表示发送信号的频域表示，$H(f)$表示系统的频率响应【对于OFDM系统，系统的频率响应是是单位脉冲响应的h(n)的DTFT】。为简化描述，下面的$Y(f)$，$X(f)$和$H(f)$分别表示为$Y$，$X$，$H$。

MMSE信道估计的目标是最小化下面的代价函数：
$$\begin{array}{rl}\mathrm{J}({\hat{\mathrm{H}}}_{\mathrm{MMSE}})& =\mathrm{E}\left\{\Vert \mathrm{e}{\Vert }^2\right\}\\ & =\mathrm{E}\left\{\Vert \mathrm{H}-{\hat{\mathrm{H}}}_{\mathrm{MMSE}}{\Vert }^2\right\}\\ & =\mathrm{E}\left\{{\Vert \mathrm{H}-\mathrm{W}{\hat{\mathrm{H}}}_{\mathrm{LS}}\Vert }^2\right\}\end{array}$$
其中，$\mathrm{H}$表示信道的真实值，${\hat{\mathrm{H}}}_{\mathrm{MMSE}}$表示MMSE信道估计值

需要特别说明的是：MMSE作为一个信道估计的方法，怎么能够在估计的时候使用真实的信道参量值呢？在实际应用过程中，由于无法获得真实信道值，因此MMSE方法无法直接应用，而是通过在仿真中预先计算加权矩阵 $W$，然后结合LS估计结果来进行MMSE信道估计
干扰管理学习日志4-------信道估计方法 LS(最小二乘)、MMSE(最小均方误差)

损失函数最小的条件为：$\mathrm{e}=\mathrm{H}-{\hat{\mathrm{H}}}_{\mathrm{MMSE}}=\mathrm{H}-\mathrm{W}{\hat{\mathrm{H}}}_{\mathrm{LS}}$ 与 ${\hat{\mathrm{H}}}_{\mathrm{LS}}$ 正交

经过一系列的计算，${\hat{H}}_{\mathrm{MMSE}}$可以表示为：

$$\begin{array}{rl}{\hat{H}}_{\mathrm{MMSE}}& =W\tilde{H}=R_{H\tilde{H}}{R}_{\tilde{H}\tilde{H}}^{-1}\tilde{H}\\ & =R_{H\tilde{H}}{\left(R_{HH}+\frac{{\sigma }_z^2}{{\sigma }_x^2}I\right)}^{-1}\tilde{H}\end{array}$$
其中， $\tilde{H}$表示LS信道估计${\hat{H}}_{\mathrm{LS}}$，$R_{H\tilde{H}}$ 是频域上真实信道$H$和 $\tilde{H}$之间的互相关矩阵,$R_{\tilde{H}\tilde{H}}$ 为 $\tilde{H}$ 的自相关矩阵，$R_{HH}$ 表示真实信道$H$的自相关矩阵，${\sigma }_z^2$ 为噪声方差, ${\sigma }_{\mathrm{x}}^2$ 为发送信号矢量方差, $I$表示单位矩阵

---

3️⃣ LMMSE(线性最小均方误差)信道估计

LMMSE信道估计是MMSE估计的一种简化形式，假设信道是线性的，信道估计的结果为：
$${\hat{H}}_{\mathrm{LMMSE}}=W\tilde{H}=R_{HH}{\left(R_{HH}+\frac{{\sigma }_z^2}{{\sigma }_x^2}I\right)}^{-1}\tilde{H}$$

其中，$R_{HH}$ 是信道的自相关矩阵，描述了信道不同时间（或频率）之间的相关性。 $\frac{{\sigma }_x^2}{{\sigma }_x^2}$ 是信噪比（SNR）的倒数。$I$ 是单位矩阵。

对于$R_{HH}$，其中的每个元素可以表示为$R_{HH}(i,j)$，代表信道在第$i$个和第$j$个子载波上的相关性。这个相关性是基于 信道的时间相关性 和 频率相关性 计算出来的：
$$R_{HH}(i,j)=\frac{1-e^{-t_{\max }(\frac{1}{t_{\mathrm{rms}}}+\mathrm{i}\cdot \frac{2\pi (i-1)}{N}})}{t_{\mathrm{rms}}\cdot \left(1-e^{-\frac{t_{\max }}{t_{\mathrm{rms}}}}\right)\cdot \left(\frac{1}{t_{\mathrm{rms}}}+\mathrm{i}\cdot \frac{2\pi (i-1)}{N}\right)}$$
$t_{\text{max }}$ 是信道延迟的最大时间
$t_{\mathrm{rms}}$ 是信道的时间相关性常数，表示信道变化的速度
$N$ 是子载波的数量
$i$ 和 $j$ 分别表示自相关矩阵中的行和列索引，对应子载波的频率索引

---

4️⃣ 代码

这里基于OFDM系统分析LS和LMMSE的代码，代码无误，可以运行：

%%%%%%%%%%%%%%%%%%%%%%LS/LMMSE信道估计OFDM仿真系统%%%%%%%%%%%%%%%%%%%
clear all；
fprintf('LS/LMMSE信道估计仿真\n');%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%
%调制方式-QPSK
%每组中2个比特看作是一个4进制的码元
modulation_mode = 4;
%导频符号与数据符号之间的间隔
pilot_interv = 8;
%OFDM符号中加入8个循环前缀
Guard = 8;
%每个OFDM符号中子载波的数量，即分组的数量
num_carrier = 64;
%OFDM符号数
num_symbol = 50;
%信噪比范围从1到20
SNR = 1:1:20;
%进行20次仿真
num_epoch = 20;
%初始化记录每个SNR和每次仿真中LS和LMMSE误码数的数组
num_bit_err_ls = zeros(length(SNR),num_epoch);
num_bit_err_lmmse = zeros(length(SNR),num_epoch);%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%
%-----------信号生成与调制-----------
for snr_index = 1:length(SNR)fprintf('信噪比为%f\n\n',SNR(snr_index));for epoch = 1:num_epoch%随机生成比特流，4进制码元，从0-3中随机生成%大小为num_carrier*num_symbolBit_Tx = randi([0,3],1,num_carrier*num_symbol);%将生成的比特流进行QPSK调制modulated_Bit = qammod(Bit_Tx,modulation_mode);%-----------导频插入与OFDM符号生成-----------%导频在插入前也需要调制，随机生成导频符号并进行QPSK调制。导频符号用于信道估计modulated_pilot = qammod(randi([0,3],num_carrier,1),modulation_mode);%计算在num_symbol个OFDM符号中，大约有多少个导频符号（ceil确保向上取整）%+1: 额外加1个导频符号，用于在数据的最后再加一个导频点，确保信道估计能够覆盖整个频域范围num_pilot = ceil(num_symbol/pilot_interv)+1;%计算OFDM帧中总的符号数（导频符号 + 数据符号）num_data = num_pilot + num_symbol;%创建一个长度为num_pilot的零向量，用于存储每个导频符号的位置pilot_index = zeros(1,num_pilot);%创建一个长度为num_symbol的零向量，用于存储数据符号的位置bit_index = zeros(1,num_symbol);%这里的num_pilot-1是因为最后一个导频位置是特殊处理的for i = 1:num_pilot-1%i-1: 这里的i是从1开始的索引，因此要先-1，让索引从0开始计算%pilot_interv+1: 每个导频之间的间隔是pilot_interv个数据符号，但是导频位置还需要向右偏移1个位置，所以是pilot_interv+1%+1: 使第一个导频符号的位置为1%pilot_interv = 8%i=1: pilot_index(1)=(1-1)*(8+1) + 1 = 1%i=2: pilot_index(2)=(2-1)*(8+1) + 1 = 10%i=3: pilot_index(3)=(3-1)*(8+1) + 1 = 19%i=4: pilot_index(4)=(4-1)*(8+1) + 1 = 28pilot_index(i) = (i-1)*(pilot_interv+1)+1;end%最后一个导频点固定在最后一个符号的位置，确保信道估计能够覆盖整个符号序列pilot_index(num_pilot) = num_data;%计算数据符号的位置%i=1: bit_index(1) = 1 + ceil(1/8) = 2%i=2: bit_index(2) = 2 + ceil(2/8) = 3%i=3: bit_index(3) = 3 + ceil(3/8) = 4%……%i=8: bit_index(8) = 8 + ceil(8/8) = 9%i=9: bit_index(9) = 9 + ceil(9/8) = 11 （因为第十个位置是导频）for i = 1:num_symbolbit_index(i) = i+ceil(i/(pilot_interv));end%插入导频data_Tx = zeros(num_carrier,num_data);%pilot_index：预先计算好的导频符号位置索引%modulated_pilot：之前已生成的导频符号（QPSK调制后的值）%repmat(modulated_pilot,1,num_pilot)：%   modulated_pilot的原始大小是 (num_carrier × 1)，即每个子载波一个导频符号%   repmat(..., 1, num_pilot)：对所有子载波上的导频符号进行重复，复制 num_pilot 次，使其适应 data_Tx 的列索引维度data_Tx(:,pilot_index) = repmat(modulated_pilot,1,num_pilot);%插入数据信号（Bit）%modulated_Bit：之前已生成的数据符号（QPSK调制后的值）%reshape(modulated_Bit, num_carrier, num_symbol)：%   modulated_Bit是一个一维向量，大小为(1 × num_carrier × num_symbol)%   reshape(..., num_carrier, num_symbol)：将其重新排列为二维矩阵，每一列对应一个完整的OFDM符号（num_carrier 个子载波）data_Tx(:,bit_index) = reshape(modulated_Bit,num_carrier,num_symbol);%-----------IFFT-----------%IFFT：频域到时域转换%ifft(data_Tx)%   data_Tx 是频域信号，大小为 (num_carrier × num_data)，表示在每个子载波上加载的符号信息（包括数据符号和导频）%   ifft(...) 进行 逆离散傅立叶变换（IDFT），将 OFDM 频域信号转换为时域信号%   由于IFFT输出的时域信号功率通常会比频域信号小（因为能量均摊到多个子载波上），所以乘以 sqrt(num_carrier) 进行功率归一化，保证时域信号的平均功率不变time_signal = sqrt(num_carrier)*ifft(data_Tx);%-----------循环前缀-----------%time_signal(num_carrier-Guard+1:num_carrier,:)：%   取OFDM符号的最后Guard个样本作为循环前缀（从 num_carrier - Guard + 1 到 num_carrier）%[time_signal(num_carrier-Guard+1:num_carrier,:); time_signal]：%    MATLAB 中的矩阵拼接操作，表示将选取的循环前缀部分和完整的OFDM符号拼接在一起，形成带有循环前缀的信号%add_cyclic_signal形状为(num_carrier + Guard) × num_data%   每一行是一个子载波的时域样本%   每一列是一个完整的OFDM符号add_cyclic_signal = [time_signal(num_carrier-Guard+1:num_carrier,:);time_signal];%-----------并串变换-----------%OFDM发送端是并行处理的%Tx_data_trans大小为：【1，(num_carrier + Guard) × num_data】Tx_data_trans = reshape(add_cyclic_signal,1,(num_carrier+Guard)*num_data);%-----------awgn信道-----------%将信号 Tx_data_trans通过加性高斯白噪声信道（AWGN） 进行传输，模拟无线信道的噪声影响%SNR(snr_index)：当前信噪比（Signal-to-Noise Ratio, SNR），即信号功率与噪声功率的比值，影响接收信号质量%'measured'：表明 MATLAB 会根据 Tx_data_trans 信号的功率 自动调整噪声功率，使得信噪比准确Rx_data = awgn(Tx_data_trans,SNR(snr_index),'measured');%-----------串并变换-----------%将串行信号恢复为(num_carrier + Guard) × num_dataRx_data_trans = reshape(Rx_data,num_carrier+Guard,num_data);%-----------去除循环前缀-----------%由于多径效应，信号在接收端可能会受到前一个符号的干扰（ISI），循环前缀的作用是消除多径效应，确保FFT变换时不会受到符号间干扰%去除循环前缀后的大小为num_carrier × num_dataRx_signal = Rx_data_trans(Guard+1:end,:);%-----------FFT-----------%fft(Rx_signal)：执行快速傅里叶变换（FFT），将时域信号转换回频域信号%/sqrt(num_carrier)：归一化 FFT 变换的输出，以确保功率与发送端匹配%data_Rx大小为num_carrier × num_datadata_Rx = fft(Rx_signal)/sqrt(num_carrier);%-----------提取导频,提取bit流-----------%modulated_pilot_Rx = data_Rx(:, pilot_index);%   提取 接收端的导频符号，用于后续信道估计%   data_Rx 是 (num_carrier × num_data) 矩阵，每列是一个 OFDM 符号，每行是不同子载波的数据%   data_Rx(:, pilot_index) 选取导频索引对应的列，得到导频信号矩阵modulated_pilot_Rx = data_Rx(:,pilot_index);%   提取接收端的数据符号，用于后续解调和误码率计算modulated_Bit_Rx = data_Rx(:,bit_index);%-----------LS信道估计-----------%modulated_pilot_Rx：接收端的导频符号%modulated_pilot：发送端的导频符号%repmat(modulated_pilot, 1, num_pilot)：将发送端的导频符号modulated_pilot重复num_pilot次，使其与modulated_pilot_Rx的维度匹配%H_ls：LS信道估计值，它是通过接收到的导频符号与发送端已知的导频符号的比值计算得到的H_ls = modulated_pilot_Rx./repmat(modulated_pilot,1,num_pilot);%进一步根据导频的估计值，利用插值法进行信道估计%pilot_index：导频符号的位置%H_ls(i, 1:num_pilot)：提取出第i个子载波上的导频符号的LS信道估计值%interp1(pilot_index, H_ls(i, 1:num_pilot), 1:num_data, 'spline')：使用spline插值方法根据导频符号的位置，对信道进行插值，估计出完整数据符号的信道响应%channel_H_ls(i, :)：得到第i个子载波的信道估计(OFDM帧中总的符号数=数据符号+导频)for i = 1:num_carrierchannel_H_ls(i,:)=interp1(pilot_index,H_ls(i,1:num_pilot),1:num_data,'spline');end%bit_index：数据符号的位置%channel_H_ls(:, bit_index)：提取数据符号对应的信道增益channel_bit_H_ls = channel_H_ls(:,bit_index);%计算通过"LS信道估计"得到的原始比特流x%modulated_Bit_Rx：接收到的数据符号%conj(channel_bit_H_ls)：信道估计的共轭复数%abs(channel_bit_H_ls).^2：信道增益的平方，用于归一化Rx_bit_esti_ls = modulated_Bit_Rx.*conj(channel_bit_H_ls)./(abs(channel_bit_H_ls).^2);%-----------LMMSE信道估计-----------%QPSK调制的参数，通常表示信号的功率beta = 1;%指系统的 时间相关性常数，用于描述信道的时间变化特性。在LMMSE中，常常使用这个值来表示信道的时间延迟或者衰减特性trms = 1e-6;%信道延迟的最大时间，这个值用于计算信道相关矩阵t_max = 8;%子载波数量N = num_carrier;%初始化一个N×N 的零矩阵，后面用来存储信道相关矩阵Rhh = zeros(N,N);%计算 信道的相关矩阵for i = 1:Nfor j = 1:NRhh(i,j) = (1-exp((-1)*t_max*((1/trms)+1i*2*pi*(i-1)/N)))./(trms*(1-exp((-1)*t_max/trms))*((1/trms)+1i*2*pi*(i-1)/N));endend%计算权重向量WW = Rhh*(eye(N)/(Rhh+(beta/SNR(snr_index)*eye(N))));%LMMSE估计中的权重矩阵H_lmmse = W*H_ls;%计算通过LMMSE信道估计得到的bit流估计值for i = 1:num_carrierchannel_H_lmmse(i,:)=interp1(pilot_index,H_lmmse(i,1:num_pilot),1:num_data,'linear');end%数据符号对应的信道增益channel_bit_H_lmmse = channel_H_lmmse(:,bit_index);%计算通过"LMMSE信道估计"得到的原始比特流Rx_bit_esti_lmmse = modulated_Bit_Rx.*conj(channel_bit_H_lmmse)./(abs(channel_bit_H_lmmse).^2);%-----------解调-----------%qamdemod：该函数用于对接收到的符号进行 QAM 解调%Rx_bit_esti_ls和Rx_bit_esti_lmmse 分别是通过LS (最小二乘法)和LMMSE(线性最小均方误差)信道估计得到的符号估计结果%modulation_mode表示所使用的调制方式%Rx_bit_esti_ls(:).'和Rx_bit_esti_lmmse(:).'是将接收的符号向量进行转置处理（从列向量变为行向量），这样可以满足qamdemod的输入要求demod_ls = qamdemod(Rx_bit_esti_ls(:).',modulation_mode);demod_lmmse = qamdemod(Rx_bit_esti_lmmse(:).',modulation_mode);%误码率计算%demod_ls(i)是通过LS信道估计解调后的第i个比特，Bit_Tx(i)是发送端第i个比特%如果解调后的比特和发送的比特不一致（即发生了误码），就将 LS 估计的误码数 (num_bit_err_ls) 增加 1for i = 1:length(Bit_Tx)if demod_ls(i)~=Bit_Tx(i)num_bit_err_ls(snr_index,epoch) = num_bit_err_ls(snr_index,epoch)+1;endif demod_lmmse(i)~=Bit_Tx(i)num_bit_err_lmmse(snr_index,epoch) = num_bit_err_lmmse(snr_index,epoch)+1;endendend
end
%%%%%%%%%%%%%%%%%%%%%%作图分析%%%%%%%%%%%%%%%%%%%
BER_ls = mean(num_bit_err_ls.')/length(Bit_Tx);
BER_lmmse = mean(num_bit_err_lmmse.')/length(Bit_Tx);figure
plot(SNR,BER_ls,'-r+',SNR,BER_lmmse,'-b*');
legend('LS','LMMSE');
title('OFDM系统LS信道估计和LMMSE信道估计比较');
xlabel('SNR');
ylabel('BER');

结果分析：

BER-SNR曲线趋势：

• LS估计：在高SNR时性能接近LMMSE，但在低SNR时误码率（BER）显著高于LMMSE。

• LMMSE估计：在低SNR时通过噪声抑制和信道统计信息利用，BER性能明显优于LS；随着SNR升高，优势逐渐减小。

性能差异原因：

• LS估计：仅利用导频点直接除法估计信道，未考虑噪声和信道相关性，对噪声敏感。

• LMMSE估计：通过信道相关矩阵$R_{HH}$ 和SNR自适应加权，抑制了噪声影响