MIMO检测之ZF,MMSE,ML算法matlab代码

最新推荐文章于 2024-07-26 17:06:46 发布

成电卷王

最新推荐文章于 2024-07-26 17:06:46 发布

阅读量5k

点赞数 6

文章标签： matlab 算法

本文链接：https://blog.csdn.net/weixin_45199852/article/details/125836647

版权

ZF和MMSE算法：

clc;
clear;
close all;

N=1000;
frame=10;
EbN0=(0:1:40);

N0=zeros(1,length(EbN0));

H=zeros(4,4,N);
W_ZF=zeros(4,4,N);
W_MMSE=zeros(4,4,N);
xzf_hat=zeros(4,1,N);
xmmse_hat=zeros(4,1,N);

y=zeros(4,1,N);

x1zf_hat=zeros(N,1);
x2zf_hat=zeros(N,1);
x3zf_hat=zeros(N,1);
x4zf_hat=zeros(N,1);

x1zfhat=zeros(N,1);
x2zfhat=zeros(N,1);
x3zfhat=zeros(N,1);
x4zfhat=zeros(N,1);

x1mmse_hat=zeros(N,1);
x2mmse_hat=zeros(N,1);
x3mmse_hat=zeros(N,1);
x4mmse_hat=zeros(N,1);

x1mmsehat=zeros(N,1);
x2mmsehat=zeros(N,1);
x3mmsehat=zeros(N,1);
x4mmsehat=zeros(N,1);

num1ZF(41)=0;
num2ZF(41)=0;
BER_ZF(41)=0;
num1MMSE(41)=0;
num2MMSE(41)=0;
BER_MMSE(41)=0;

h11=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h12=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h13=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h14=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h21=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h22=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h23=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h24=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h31=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h32=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h33=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h34=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h41=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h42=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h43=(randn(1,N)+1j*randn(1,N))/sqrt(2);
h44=(randn(1,N)+1j*randn(1,N))/sqrt(2);

x1=round(rand(N,1));
x2=round(rand(N,1));
x3=round(rand(N,1));
x4=round(rand(N,1));

u1=2*x1-1;
u2=2*x2-1;
u3=2*x3-1;
u4=2*x4-1;

for t=1:N
    H(:,:,t)=[h11(t),h12(t),h13(t),h14(t);h21(t),h22(t),h23(t),h24(t);h31(t),h32(t),h33(t),h34(t);h41(t),h42(t),h43(t),h44(t)];
end

for i=1:length(EbN0)
    berZF=0;
    berMMSE=0;
    berML=0;
    N0(i)=1/(10^((EbN0(i))/10));
    for m=1:frame
        z1=sqrt(N0(i)/(2))*(randn(size(u1))+1j*randn(size(u1)));
        z2=sqrt(N0(i)/(2))*(randn(size(u2))+1j*randn(size(u2)));
        z3=sqrt(N0(i)/(2))*(randn(size(u3))+1j*randn(size(u3)));
        z4=sqrt(N0(i)/(2))*(randn(size(u4))+1j*randn(size(u4)));
        
        y1 = u1.*(h11).'+u2.*(h12).'+u3.*(h13).'+u4.*(h14).'+z1.';
        y2 = u1.*(h21).'+u2.*(h22).'+u3.*(h23).'+u4.*(h24).'+z2.';
        y3 = u1.*(h31).'+u2.*(h32).'+u3.*(h33).'+u4.*(h34).'+z3.';
        y4 = u1.*(h41).'+u2.*(h42).'+u3.*(h43).'+u4.*(h44).'+z4.';
        
        for q=1:N
            y(:,:,q)=[y1(q);y2(q);y3(q);y4(q)];
        end
        %ZF================================================================
        for n=1:N
            W_ZF(:,:,n)= inv(conj(H(:,:,n).')*H(:,:,n))*(conj(H(:,:,n).'));
        end
        
        for p=1:N
            xzf_hat(:,:,p)=W_ZF(:,:,p)*y(:,:,p);
        end
        
        for o=1:N
            x1zf_hat(o)=xzf_hat(1,1,o);
            x2zf_hat(o)=xzf_hat(2,1,o);
            x3zf_hat(o)=xzf_hat(3,1,o);
            x4zf_hat(o)=xzf_hat(4,1,o);
        end
        
        for k=1:N
            if x1zf_hat(k) < 0
                x1zfhat(k) = 0;
            else
                x1zfhat(k) = 1;
            end
            
            if x2zf_hat(k) < 0
                x2zfhat(k) = 0;
            else
                x2zfhat(k) = 1;
            end
            
            if x3zf_hat(k) < 0
                x3zfhat(k) = 0;
            else
                x3zfhat(k) = 1;
            end
            
            if x4zf_hat(k) < 0
                x4zfhat(k) = 0;
            else
                x4zfhat(k) = 1;
            end
        end
        
        [num1ZF, rat1ZF] = biterr(x1zfhat, x1);
        [num2ZF, rat2ZF] = biterr(x2zfhat, x2);
        [num3ZF, rat3ZF] = biterr(x3zfhat, x3);
        [num4ZF, rat4ZF] = biterr(x4zfhat, x4);
        berZF=berZF+(rat1ZF+rat2ZF+rat3ZF+rat4ZF)/4;
        
        %MMSE==============================================================
        for n=1:N
            W_MMSE(:,:,n)= inv(conj(H(:,:,n).')*H(:,:,n)+N0(i)*eye(4))*(conj(H(:,:,n).'));
        end
        
        for p=1:N
            xmmse_hat(:,:,p)=W_MMSE(:,:,p)*y(:,:,p);
        end
        
        for o=1:N
            x1mmse_hat(o)=xmmse_hat(1,1,o);
            x2mmse_hat(o)=xmmse_hat(2,1,o);
            x3mmse_hat(o)=xmmse_hat(3,1,o);
            x4mmse_hat(o)=xmmse_hat(4,1,o);
        end
        
        for k=1:N
            if x1mmse_hat(k) < 0
                x1mmsehat(k) = 0;
            else
                x1mmsehat(k) = 1;
            end
            if x2mmse_hat(k) < 0
                x2mmsehat(k) = 0;
            else
                x2mmsehat(k) = 1;
            end
            if x3mmse_hat(k) < 0
                x3mmsehat(k) = 0;
            else
                x3mmsehat(k) = 1;
            end
            if x4mmse_hat(k) < 0
                x4mmsehat(k) = 0;
            else
                x4mmsehat(k) = 1;
            end
        end
        [num1MMSE, rat1MMSE] = biterr(x1mmsehat, x1);
        [num2MMSE, rat2MMSE] = biterr(x2mmsehat, x2);
        [num3MMSE, rat3MMSE] = biterr(x3mmsehat, x3);
        [num4MMSE, rat4MMSE] = biterr(x4mmsehat, x4);
        berMMSE=berMMSE+(rat1MMSE+rat2MMSE+rat3MMSE+rat4MMSE)/4;
        
    end
    BER_ZF(i)=berZF/m;
    BER_MMSE(i)=berMMSE/m;
end

figure
semilogy(EbN0,BER_ZF,'o-b',EbN0,BER_MMSE,'*-r');
title('误码率');
xlabel('EbN0(dB)');
ylabel('BER');
legend('ZF','MMSE');
grid on;

ML算法：

close all;
clear;
clc;
Nt = 2;
Nr = 2;
Len = 10000;
M= 4;
bitsPerSymbol = log2(M) ;
bitsTotal = bitsPerSymbol * Nt * Len ;
ALP = [-1+1i  ;-1-1* 1i ;1+1i ;1-1i] * sqrt(3/(2*(M-1)))  ; % QPSK符号
SNRVect = 0:2:20 ;
SER = zeros( 1 ,length(SNRVect) )  ;
for snrLoop =  1: length (SNRVect)
    fprintf(' snrLoop = %d \n',snrLoop) ;
    % 计算信噪比
    SNRdB = SNRVect(snrLoop) ;
    SNR  = 10^(SNRdB/10) ;
    % 此处我们假设每根天线发送的信号能量为1   
    % N0 表示噪声的能量
    N0 = Nt * 1 / SNR ;   
    % 产生随机比特流
    send = round(rand(1,bitsTotal)) ;
    sendReshape = reshape(send , Nt * bitsPerSymbol , Len ) ;
    symbols = zeros(Nt, Len) ;
    symbolsModulation = zeros(Nt ,Len) ;
    % 发送端QPSK调制
    fprintf('\n正在调制并发送 \n') ;
    for sendLoop = 1 :Len
        for antennaLoop = 1: Nt
            twoBits = sendReshape ( 2* antennaLoop-1 :2* antennaLoop, sendLoop);
            symbols(antennaLoop ,sendLoop) = 2 * twoBits(1) + 1 * twoBits(2) +1;
            symbolsModulation(antennaLoop ,sendLoop) = ALP(symbols(antennaLoop ,sendLoop));
        end
    end
    % 经过信道
    H = (randn(Nr, Nt) + sqrt(-1) * randn(Nr , Nt))/ sqrt(2)  ;            % 除以sqrt(2) 保证信道的能量为1
    n = (randn(Nr , Len) + sqrt(-1)*randn(Nr, Len )) * sqrt(N0) /sqrt(2) ; % 保证噪声 n ~N(0,N0)
    y = H* symbolsModulation + n ;
    % 接收端
    fprintf('\n正在译码 \n');
    symbolsDetect = zeros(2 , Len);
    for sendLoop = 1 :Len
         fprintf('%d %d \n',snrLoop ,sendLoop);
        % 得到M*M种组合
        dml = zeros(M ,M ,Len);
        for i1 = 1:M
            for i2 = 1:M
                dml(i1, i2, sendLoop) = (norm(y(:,sendLoop) - H*[ALP(i1);ALP(i2)])).^2; % 共16种组合
            end
        end
        dml_min = min(reshape(dml(:,:,sendLoop).',1,M*M));
        % 找到最合适的然后译码
        for j1 = 1:M 
            for j2 = 1:M
                if dml(j1 , j2, sendLoop) == dml_min 
                    symbolsDetect(1, sendLoop) = j1;
                    symbolsDetect(2, sendLoop) = j2;
                end
            end
        end

    end
    % 一个信噪比结束  统计SER
    SER(snrLoop) = sum(reshape(symbolsDetect,1,Len*Nt) ~= reshape(symbols,1,Len*Nt))/(Len*Nt);
end
figure
semilogy(0:2:20,SER);