1. MMSE信道估计
在线性最小均方误差估计（LMMSE）中，采用的准则函数是统计均方误差，首先将待定信道估计值表示为：
其中W是待定的加权系数矩阵，这样就是要估计这些权系数，LMMSE估计的准则函数如下：
对W求偏导，并令其等于0，利用正交性原理，可得：
即
从而得到G的LMMSE为
最后得到H的估计为
下图是在EVA信道条件下，信道估计技术与均衡技术的性能比较图，其中调制方式为16-QAM，编码采用1/2码率的Turbo码。
由图可见，在多径信道条件下，LMMSE信道估计性能还是明显要优于LS信道估计。也可以看到在低信噪比条件下，性能并没有很大改善，这时噪声的影响占了主要因素。如果在AWGN信道条件下，由于Turbo码的存在，LS信道估计的性能在信噪比12dB的情况下已经不出现误码率，此时两种信道估计的曲线相差不大，因此在本系统的信道条件下（信道条件为室内信道，信噪比14dB以上），采用LS信道估计已足够。另外，这里由于天线间的干扰不存在，所以MMSE均衡和ZF均衡性能曲线几乎重合。
考虑到可能存在的定时误差，将导致相位偏移，相位偏差是关于k的线性函数。采用跟踪导频做信道估计时，忽略相位噪声的影响，采用线性内插时可将所有的相位偏移准确的估计出来，此时在频域采用线性内插最合适。在室内信道条件下，能保证信噪比大于14dB，此时信道编码可以完全正确译码，因此频域的信道估计方法采用LS算法。
在白高斯信道条件下，如果不存在定时误差，估计出的信道应该是信道响应的实部近似平坦，虚部为零；若存在定时误差，会使得估计出的信道有一定的相位偏转，两种情况下信道估计结果如下图所示
上图为信道估计值的星座图，下图为时域图，时域图中上图为频域内插之后的信道估计值的实部，下图为虚部，前600个子载波为正频率，后600个子载波为负频率。第一幅图为定时准确的情况，可见信道估计值仍然存在相位噪声；第二幅图是定时点出现误差的情况，此时引入的相位偏差会叠加到相位噪声上，如果定时偏差在3个样点以内，将引入小于2pi的相位偏转，该图的定时偏差刚好为3个样点。
实际信道为室内信道，时域上信道变化缓慢，因此可采用最邻近插值的方法估计其他未插导频的符号。此时认为在一个时隙内的信道响应近似不变。
因此信道估计算法频域采用线性插值，时域采用邻近插值。位于天线端口2上的导频处未插入任何导频信息，可以利用此特性来做信噪比估计，即使用每个时隙第6个符号上的导频计算信号功率，第7个符号上的导频计算噪声功率。
考虑到存在一定的相位噪声，使得每个符号的星座图有一定的相位偏转，一个主要的相位噪声是残余频偏的影响。假设在频偏补偿之后残余频偏还有f，则使用上述时域信道内插方法使得这些频偏会累计到一个时隙上。假设最后的残余频偏为f_delta，则一个时隙内每个符号残余相偏为
由于时域上信道估计是直接扩展得到的，所以残余相偏会因为时间的累计而增加。累计的间隔为一个OFDM符号长度（加CP）。整个下行子帧所有OFDM符号的相位偏转如下图所示：
每一个OFDM符号经历相同的相位偏转，当频偏补偿效果理想残余频偏为10Hz时，最大相偏约为0.01pi，实际信道下一个时隙的每个符号的星座图相偏如下所示
2. 请问如何做信道估计中的MMSE估计
MMSE估计就是最小均方误差估计，通过求得一个合适的信道冲击响应（CIR），使得通过CIR计算出的接收数据与实际数据的误差的均方和最小。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
我上个月刚做过基于块状导频信息的LTE物理层上行信道的频域信道估计以及信道均衡。
部分算法如下（以下是基于单载波的）
假设循环前缀已经消除了实践弥散信道带来的符号间干扰，保证了子载波之间的正交性。并且信道为慢衰落信道，在一个OFDM符号内，可以认为保持不变。
均衡器接收到的信号可以表示为
y(t)=x(t)*h(t)+n(t)
y(t)为均衡器接收到的信号，h(t)为系统等效的冲击响应，x(t)为原始的输入信号，n(t)为系统中的噪声。
信道估计的任务就是在已知发送参考信息的情况下，对接受到的参考信息进行分析，选择合适的算法得到参考信息的信道冲击响应，即h(t),而数据信息的信道冲击响应则可以通过插值得到。
1) 最小二乘估计（LS）
该算法的目的是
有正交性原理，则可得LS估计
该估计为无偏估计，每估计一个新到衰落系数只需一次乘法，缺点是受噪声影响较大。
2) 线性最小均方误差估计（MMSE）
LMMSE估计属于统计估计，需要对信道的二阶统计量进行估计，利用信道相关性可以置信道噪声提高估计性能。以最小均方误差（MMSE）为准则，如下式：
为了降低计算的复杂度，一般将 用它的期望值 代替，信道性能不会产生明显恶化，则上式可变为
其中 为一个仅与调试的星座的大小有关的值， 为平均信噪比。
该算法的复杂度较高，随着X的改变， 须不断更新。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
不知道你的是物理模型和数据结构是什么样的，频域估计还是时域估计，基于导频信息还是盲信道估计？
____________________________________________
有点悲剧，Word里面的公式我不知道怎么插进来
3. MATLAB 程序详解（关于波束形成）
你这里有两个程序，第二个程序与第一个实质上是一样的，区别就是信号与导向矢量的写法有点不同，这里我就不注释了。还有，我下面附了一段我自己的写的程序，里面有SIM算法。G-S正交化算法等。是基于圆阵形式的，你的算法是基于线阵的，他们程序上的区别在于导向矢量的不同。我的算法是某项目中的，保证好使。建议学习波束形成技术，注意把程序分块，例如分成，求导向矢量；最优权值；形成波束等等。
程序如下：
4单元均匀线阵自适应波束形成图
clear
clc
format long;
v=1;
M=4;
N=1000;%%%%%%%快拍数
f0=21*10^3;%%%%%%%%%%%信号与干扰的频率
f1=11*10^3;
f2=15*10^3;
omiga0=2*pi*f0;%%%%%%%信号与干扰的角频率
omiga1=2*pi*f1;
omiga2=2*pi*f2;
sita0=0.8; %信号方向
sita1=0.4; %干扰方向1
sita2=2.1; %干扰方向2
for t=1:N %%%%%%%%%%%%信号
adt(t)=sin(omiga0*t/(N*f0));
a1t(t)=sin(omiga1*t/(N*f1));
a2t(t)=sin(omiga2*t/(N*f2));
end
for i=1:M %%%%%%%%%%%%信号的导向矢量：线阵的形式
ad(i,1)=exp(j*(i-1)*pi*sin(sita0));
a1(i,1)=exp(j*(i-1)*pi*sin(sita1));
a2(i,1)=exp(j*(i-1)*pi*sin(sita2));
end
R=zeros(M,M);
for t=1:N
x=adt(t)*ad+a1t(t)*a1+a2t(t)*a2; %阵列对信号的完整响应
R=R+x*x';%信号的协方差矩阵
end
R=R/N;%%%%%%%%%协方差矩阵，所有快拍数的平均
miu=1/(ad'*inv(R)*ad);%%%%%%这个貌似是LMS算法的公式，具体我记不太清，这里是求最优权值，根据这个公式求出，然后加权
w=miu*inv(R)*ad;
%%%%%%形成波束%%%%%%%%%%%%%%%%%%%
for sita=0:pi/100:pi
for i=1:M
x_(i,1)=exp(j*(i-1)*pi*sin(sita));
end
y(1,v)=w'*x_;%%%%%%%对信号进行加权，消除干扰
v=v+1;
end
y_max=max(y(:));%%%%%%%%%%%%%%%归一化
y_1=y/y_max;
y_db=20*log(y_1);
sita=0:pi/100:pi;
plot(sita,y)
Xlabel（‘sitaa’）
Ylabel（‘天线增益db’）
4单元均匀线阵自适应波束形成
目标
clear
clc
format long;
v=1;
M=4;阵元数
N=100;
f0=21*10^3;
omiga0=2*pi*f0;
sita0=0.6;%信号方向
for t=1:N
adt(t)=sin(omiga0*t/(N*f0));
end
for i=1:M
ad(i,1)=exp(j*(i-1)*pi*sin(sita0));
end
R=zeros(4,4);
r=zeros(4,1);
for t=1:N
x=adt(t)*ad;
R=R+x*x.';
end
R=R/N;
miu=1/(ad.'*inv(R)*ad);
w=miu*inv(R)*ad;
for sita=0:pi/100:pi/2
for i=1:M
a(i,1)=exp(j*(i-1)*pi*sin(sita));
end
y(1,v)=w.'*a;
v=v+1;
end
sita=0:pi/100:pi/2;
plot(sita,y)
xlabel('sita')
ylabel('天线增益’)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%我的程序%%%%%%%%%%%%%%%
function jieshousignal
%期望信号数：1个
%干扰信号数：4个
%信噪比已知
%干燥比已知
%方位角已知
clc;
clear all;
close all;
%//参数设置===========================================
www1=0;
www2=0;
www3=0;
% for rrr=1:16000
signal_num=1; %signal number
noise_num=5; %interference number
R0=0.6; %圆的半径
SP=2000; %Sample number
N=8; %阵元数
snr=-10; %Signal-to-Noise
sir1=10; %Signal-to-Interference one
sir2=10; %Signal-to-Interference two
sir3=10; %Signal-to-Interf
sir4=10;
sir5=10;
%//================noise Power-to-signal Power====================
factor_noise_1=10.^(-sir1/10);
factor_noise_2=10.^(-sir2/10);
factor_noise_3=10.^(-sir3/10);
factor_noise_4=10.^(-sir4/10);
factor_noise_5=10.^(-sir5/10);
factor_noise_targe=10.^(-snr/10);
% //======================== ===============
d1=85*pi/180;%%干扰1的方位角
d2=100*pi/180;%干扰2的方位角
d3=147*pi/180;%干扰3的方位角
d4=200*pi/180;%干扰4的方位角
d5=250*pi/180;%干扰5的方位角
d6=150*pi/180;%目标的方位角
e1=15*pi/180;%%干扰1的俯仰角
e2=25*pi/180;%干扰2的俯仰角
e3=85*pi/180;%干扰3的俯仰角
e4=50*pi/180;%干扰4的俯仰角
e5=70*pi/180;%干扰5的俯仰角
e6=85*pi/180;%目标的俯仰角
% //====================目标信号==========================
t=1:1:SP;
fc=2e7;
Ts=1/(3e10);
S0=5*cos(2*pi*fc*t*Ts);%目标信号
for kk=1:N
phi_n(kk)=2*pi*(kk-1)/N;
end
%//====================操纵矢量==========================================
A=[conj(exp(j*2*pi*R0*cos(d6-phi_n)*sin(e6)));conj(exp(j*2*pi*R0*cos(d1-phi_n)*sin(e1)));conj(exp(j*2*pi*R0*cos(d2-phi_n)*sin(e2)));conj(exp(j*2*pi*R0*cos(d3-phi_n)*sin(e3)));conj(exp(j*2*pi*R0*cos(d4-phi_n)*sin(e4)));conj(exp(j*2*pi*R0*cos(d5-phi_n)*sin(e5)))]';
A1=[conj(exp(j*2*pi*R0*cos(d1-phi_n)*sin(e1)));conj(exp(j*2*pi*R0*cos(d2-phi_n)*sin(e2)));conj(exp(j*2*pi*R0*cos(d3-phi_n)*sin(e3)));conj(exp(j*2*pi*R0*cos(d4-phi_n)*sin(e4)));conj(exp(j*2*pi*R0*cos(d5-phi_n)*sin(e5)))]';
% //==========================================================Power of the interference
% // depending on the signal power and SIR
Ps1=0;
Ps2=0;
Ps3=0;
Ps4=0;
Ps5=0;
S1=zeros(1,SP);
S2=zeros(1,SP);
S3=zeros(1,SP);
S4=zeros(1,SP);
S5=zeros(1,SP);
Ps0=S0*S0'/SP; % signal power
Ps1=Ps0*factor_noise_1;
Ps2=Ps0*factor_noise_2;
Ps3=Ps0*factor_noise_3;
Ps4=Ps0*factor_noise_4;
Ps5=Ps0*factor_noise_5;
% //==========================干扰信号的随机包络=========================
S1=normrnd(0,sqrt(Ps1/2),1,SP)+j*normrnd(0,sqrt(Ps1/2),1,SP);
S2=normrnd(0,sqrt(Ps2/2),1,SP)+j*normrnd(0,sqrt(Ps2/2),1,SP);
S3=normrnd(0,sqrt(Ps3/2),1,SP)+j*normrnd(0,sqrt(Ps3/2),1,SP);
S4=normrnd(0,sqrt(Ps4/2),1,SP)+j*normrnd(0,sqrt(Ps4/2),1,SP);
S5=normrnd(0,sqrt(Ps5/2),1,SP)+j*normrnd(0,sqrt(Ps5/2),1,SP);
%//
S=[S0;S1;S2;S3;S4;S5];
SS1=[S1;S2;S3;S4;S5];
X=A*S;%信号加干扰
XX2=A1*SS1; %接收到的干扰
Pw_noise=sqrt(Ps0*factor_noise_targe/2);
a1=randn(N,SP);
a2=randn(N,SP);
a1=a1/norm(a1);
a2=a2/norm(a2);
W=Pw_noise*(a1+sqrt(-1)*a2);
X=X+W;
% //--------------------------SMI算法----------------------------------------
Rd=X*S0'/SP;
R=X*X'/(SP*1);
Wc_SMI=pinv(R)*Rd./(Rd'*pinv(R)*Rd);%权向量
Wc_SMI=Wc_SMI/norm(Wc_SMI);
Y_SMI=Wc_SMI'*X; %SMI算法恢复出来的信号
%//-------------------------------------GS算法------------------
m=1;
for i=1:400:2000
X2(:,m)=XX2(:,i);
m=m+1;
end
a=zeros(1,8);
phi_n=zeros(1,8);
phi=0:pi/180:2*pi;
theta=0:pi/180:pi/2;
for kk=1:8
a(kk)=1;
phi_n(kk)=2*pi*(kk-1)/8;
end
x1=zeros(1,8);
x2=zeros(1,8);
x3=zeros(1,8);
x4=zeros(1,8);
x5=zeros(1,8);
x1=X2(:,1)';
x2=X2(:,2)';
x3=X2(:,3)';
x4=X2(:,4)';
x5=X2(:,5)';
Z1=x1;
Z1_inner_proct=Z1.*conj(Z1);
Z1_mode=sqrt(sum(Z1_inner_proct));
Y1=Z1./Z1_mode;
Inner_proct=sum(x2.*conj(Y1));
Z2=x2-Inner_proct*Y1;
Z2_inner_proct=sum(Z2.*conj(Z2));
Z2_mode=sqrt(Z2_inner_proct);
Y2=Z2./Z2_mode;
Inner_proct1=sum(x3.*conj(Y1));
Inner_proct2=sum(x3.*conj(Y2));
Z3=x3-Inner_proct1*Y1-Inner_proct2*Y2;
Z3_inner_proct=sum(Z3.*conj(Z3));
Z3_mode=sqrt(Z3_inner_proct);
Y3=Z3./Z3_mode;
Inner_proct1_0=sum(x4.*conj(Y1));
Inner_proct2_0=sum(x4.*conj(Y2));
Inner_proct3_0=sum(x4.*conj(Y3));
Z4=x4-Inner_proct1_0*Y1-Inner_proct2_0*Y2-Inner_proct3_0*Y3;
Z4_inner_proct=sum(Z4.*conj(Z4));
Z4_mode=sqrt(Z4_inner_proct);
Y4=Z4./Z4_mode;
Inner_proct1_1=sum(x5.*conj(Y1));
Inner_proct2_1=sum(x5.*conj(Y2));
Inner_proct3_1=sum(x5.*conj(Y3));
Inner_proct4_1=sum(x5.*conj(Y4));
Z5=x5-Inner_proct1_1*Y1-Inner_proct2_1*Y2-Inner_proct3_1*Y3-Inner_proct4_1*Y4;
Z5_inner_proct=sum(Z5.*conj(Z5));
Z5_mode=sqrt(Z5_inner_proct);
Y5=Z5./Z5_mode;
%Y1
%Y2
%Y3
%Y4
%Y5
w0=zeros(1,8);
w=zeros(1,8);
for mm=1:8;
w0(mm)=exp(-j*2*pi*R0*cos(d6-phi_n(mm))*sin(e6));
end
dd1=sum(w0.*conj(Y1))*Y1;
dd2=sum(w0.*conj(Y2))*Y2;
dd3=sum(w0.*conj(Y3))*Y3;
dd4=sum(w0.*conj(Y4))*Y4;
dd5=sum(w0.*conj(Y5))*Y5;
w=w0-dd1-dd2-dd3-dd4-dd5;
Wc_GS=w;
Wc_GS=Wc_GS/(norm(Wc_GS));
Y_GS=Wc_GS*X; %GS算法恢复出来的图像
%//----------------------------------MMSE算法-----------------------
Rd=X*S0'/SP;
R=X*X'/(SP*1);
Wc_MMSE=pinv(R)*Rd;
Wc_MMSE=Wc_MMSE/norm(Wc_MMSE);
Y_MMSE=Wc_MMSE'*X; %MMSE算法恢复出来的信号
S0=S0/norm(S0);
Y_GS=Y_GS/norm(Y_GS);
Y_SMI=Y_SMI/norm(Y_SMI);
Y_MMSE=Y_MMSE/norm(Y_MMSE);
% figure(1)
% plot(real(S0));
% title('原始信号');
% xlabel('采样快拍数');
% ylabel('信号幅度');
% figure(2)
% plot(real(Y_SMI));
% title('运用SMI算法处理出的信号');
% xlabel('采样快拍数');
% ylabel('信号幅度');
% figure(3)
% plot(real(Y_GS));
% title('运用G-S算法处理出的信号');
% xlabel('采样快拍数');
% ylabel('信号幅度');
% figure(4)
% plot(real(Y_MMSE));
% for i=1:SP
% ss(i)=abs(S0(i)-Y_SMI(i))^2;
% end
% q_1=mean(ss);
% for i=1:SP
% ss1(i)=abs(S0(i)-Y_GS(i))^2;
% end
% q_2=mean(ss1);
% for i=1:SP
% ss2(i)=abs(S0(i)-Y_MMSE(i))^2;
% end
% q_3=mean(ss2);
%
% www1=www1+q_1;
% www2=www2+q_2;
% www3=www3+q_3;
% end
% www1/16000
% www2/16000
% www3/16000
phi=0:pi/180:2*pi;
theta=0:pi/180:pi/2;
%
% % //------------------------ 形成波束-----------------------------------------
F_mmse=zeros(91,361);
F_smi=zeros(91,361);
F_gs=zeros(91,361);
for mm=1:91
for nn=1:361
p1=sin(theta(mm));
p2=cos(phi(nn));
p3=sin(phi(nn));
q1=sin(e6);
q2=cos(d6);
q3=sin(d6);
for hh=1:8
w1=cos(phi_n(hh));
w2=sin(phi_n(hh));
zz1=q2*w1+q3*w2;
zz2=p2*w1+p3*w2;
zz=zz2*p1-zz1*q1;
F_mmse(mm,nn)= F_mmse(mm,nn)+conj(Wc_MMSE(hh))*(exp(j*2*pi*R0*(zz2*p1)));
F_smi(mm,nn)=F_smi(mm,nn)+conj(Wc_SMI(hh))*(exp(j*2*pi*R0*(zz2*p1)));
F_gs(mm,nn)=F_gs(mm,nn)+conj((Wc_GS(hh))')*(exp(j*2*pi*R0*(zz2*p1)));
end
end
end
F_MMSE=abs(F_mmse);
F_SMI=abs(F_smi);
F_GS=abs(F_gs);
figure(5)
mesh(20*log10(F_MMSE))
figure(6)
mesh(20*log10(F_SMI))
title('SMI算法波束形成图');
xlabel('方位角');
ylabel('俯仰角');
zlabel('幅度/dB');
figure(7)
mesh(20*log10(F_GS))
title('G-S算法波束形成图');
xlabel('方位角');
ylabel('俯仰角');
zlabel('幅度/dB');
4. 请问下谁有MIMO系统中的ML、MMSE或者ZF等检测算法的matlab代码呀小妹急需使用，非常感谢，非常感谢呀！
pudn上面随便咦搜就一大堆，
不过看别人的代码真的很费劲，有些人的思维你真的没法理解。
还不如自己搞懂公式了自己来。