mainar1.m

此程序为基于KALMAN滤波的信道估计联合符号检测算法。 mainar1中信道模型为AR1模型。

clear;clc;

N=100; %total number of symbol pairs to be transmitted (should be at least 10 times more than expected 1/min(BER))
M=4; %PSK order (must be a power of 2): 2, 4, 8 etc'
Tx=2; %number of Tx elements, must be 2
Rx=1; %number of Rx elements
SNR=0:3:30; %SNR in dB, average received power at one Rx element over the average noise power at that element
tr=floor(N/12);
aver=1;

BER_temp=zeros(1,length(SNR));
for m=1:aver
Alpha=0.998;
delta=1-Alpha^2;
h(1,:)=(randn(1,2)+j*randn(1,2))/sqrt(2);     
p=1;
for n=2:2*N+2*tr+1;
    v2=sqrt(delta/2)*(randn(1,2)+j*randn(1,2)); 
    h(n,:)=Alpha*h(n-1,:)+v2;                 
    ch(p,:)=h(n,:); 
    p=p+1;
end   
for k=1:length(SNR);
    %Toss pairs of uniformly distributed MPSK symbols with power 1/2
    s=[];z=[];
    Ctr=[0 2].';
    train=exp(j*2*pi/M*Ctr)*2;
    C=floor(M*rand(2,N)); %transmitted alphabet
    st=exp(j*2*pi/M*C)/sqrt(2);
    %Simulate equivalent matrix of impulse noise
    %Noise power caculation
    snr=10^(SNR(k)/10); %just translate SNR from dB to times
    
    sig=0.5/snr; %the sigma square of the noise
    
    Ns=sqrt(sig)*(randn(2*Rx,N+tr)+j*randn(2*Rx,N+tr)); %noise matrix 
        
    t=1;q=1;p=0;
    for n=2:2:2*N+2*tr
        p=p+1;
        h1=h(n,:);
        h2=h(n+1,:);
        H=[h1(1) h1(2); h2(2)' -h2(1)'];
        if mod(p,13)==0
           sr(:,t)=H*train+Ns(:,t);
           sr1(:,t)=train(1)*h1(1)+train(2)*h1(2)+Ns(1,t);
           sr2(:,t)=-train(2)'*h2(1)+train(1)'*h2(2)+Ns(2,t);
           t=t+1;
           p=0;
        else
           sr(:,t)=H*st(:,q)+Ns(:,t);
           sr1(:,t)=st(1,q)*h1(1)+st(2,q)*h1(2)+Ns(1,t);
           sr2(:,t)=-st(2,q)'*h2(1)+st(1,q)'*h2(2)+Ns(2,t);
           t=t+1;
           q=q+1; 
        end    
   end
    W=2*sig;  
    V=[delta 0;0 delta];
    A=[Alpha 0;0 Alpha];
    z.h2=h(1,:).';
    z.P=[0 0;0 0]';
    p=0;x=1;t=1;
    for n=1:N+tr;
        p=p+1;
        sy=sr(:,n);
        sy1=sr1(:,n);
        sy2=sr2(:,n);
        s.P0 = A * z.P * A' + V;  
        s.h10 = Alpha*z.h2;
        z.h20=Alpha^2*z.h2;
        if mod(p,13)==0
           s.d1=train.';
           z.d2=[-train(2)' train(1)'];
        else
           H_est=[s.h10(1) s.h10(2);z.h20(2)' -z.h20(1)'];
           sd1=H_est'*sr(:,n); %received symbols
           %ML detection
           ang=angle(sd1); %received phase angles
           B1=mod(round(ang/(2*pi/M)),M);
           s.d1=(exp(j*2*pi/M*B1)/sqrt(2)).';
           z.d2=[-s.d1(2)' s.d1(1)'];
        end
        for i=1:10;
            s=kalmanf1(s,W,sy1); % perform a Kalman filter iteration
            z.h20 = A * s.h1;           
            z.P0 = A * s.P * A' + V;                              
            z=kalmanf2(z,W,sy2); 
            % refined decoder
            if mod(p,13)~=0
               H=[s.h1(1) s.h1(2); z.h2(2)' -z.h2(1)'];
               % ML decoder
               sd=H'*sy;
               ang=angle(sd); %received phase angles
               B=mod(round(ang/(2*pi/M)),M); %received alphabet
               s.d1=(exp(j*2*pi/M*B)/sqrt(2)).';
               z.d2=[-s.d1(2)' s.d1(1)'];
            end
        end   
        chann_est(:,t)=s.h1;
        chann_est(:,t+1)=z.h2;
        t=t+2;
        if mod(p,13)==0
           p=0;
        else
           D(:,x)=B;
           x=x+1;
        end
    end
    %BER estimation (for the Gray code constellations)
    BER(:,k)=sum(sum(xor(D-C,0)))/2/N/log2(M);  
end %k (SNR)
BER_temp=BER_temp+BER;
end
BER=BER_temp/aver;
figure(1)
semilogy(SNR, BER, 'g-v');
xlabel('SNR, [dB]');
ylabel('BER');
grid on
hold on