mmse2.m

是一个对mimo连续相位调制系统的均衡程序

function L_all = mmse2( noisy_sig, mtchd_fltr, D, Ns, H, L, sigma, L_a, M)
% This script used for demodulating CPM signals
% Input: 
%        noisy_sig: Received signals
%        mtchd_fltr: Matched filter
%        D: The duration of the pulse
%        Ns: Number of samples per symbol
%        M: samples rate  
%        H: Estimated channel coefficients
%        sigma : Deviation of  noise
% Output:
%        L_all: logarithm likelihood ratio per bit

mod_sym_num = length(noisy_sig)/Ns - L + 1;
noisy_sig = [noisy_sig.';zeros((D-L)*Ns,M)];
R = zeros(M, mod_sym_num);
r = zeros(M, mod_sym_num/2);
dhat = zeros(1, mod_sym_num);
Mean = zeros(1, mod_sym_num);

for k = 1 : mod_sym_num/2
    for l = 1 : M
        R(l,2*k-1) = sum(noisy_sig((2*k-2)*Ns+1:(2*k-2+D)*Ns,l).*mtchd_fltr); % Info. after Matched filter 
        R(l,2*k)   = sum(noisy_sig((2*k-1)*Ns+1:(2*k-1+D)*Ns,l).*mtchd_fltr); % Info. after Matched filter
        r1(l,k) = R(l,2*k-1);
        r2(l,k) = R(l,2*k);
    end
end

rr1 = [r1 zeros(M,1)];
rr2 = [zeros(M,1) r2];
RR1 = reshape(rr1,M*(mod_sym_num/2+1),1);
RR2 = reshape(rr2,M*(mod_sym_num/2+1),1);
% rr = [zeros(M,1) r zeros(M,1)];
% RR = reshape(rr,M*(mod_sym_num/2+2),1);

L_a = [zeros(1,2) L_a zeros(1,2)]'; % priori info.
dd  = tanh(-L_a/2);                 % estimated info.

% With MMSE FIR
for l=1:D
    c(l) = sum(mtchd_fltr(1:(D-l+1)*Ns).*mtchd_fltr((l-1)*Ns+1:D*Ns)); % Matched coefficients
end

% Sigma = 2 * sigma^2;
Sigma = sigma^2;

j = sqrt(-1);
for l = 1 : M
    h1(l,:) = [H(l,2)*c(2)  H(l,1)*c(2)   j*H(l,1)*c(1)+H(l,2)*c(2)  H(l,1)*c(2)-j*H(l,2)*c(1)]/sqrt(2*M); % ISI coefficients
    h2(l,:) = [j*H(l,1)*c(2)+H(l,2)*c(1)  H(l,1)*c(1)-j*H(l,2)*c(2)  j*H(l,1)*c(2)  -j*H(l,2)*c(2)]/sqrt(2*M); % ISI coefficients
end
HH1 = [h1 zeros(M,2);zeros(M,2) h1];
HH2 = [h2 zeros(M,2);zeros(M,2) h2];
% HH = [h zeros(M,4); zeros(M,2) h zeros(M,2); zeros(M,4) h];

o_3 = [zeros(1,2) 1 zeros(1,3)]'; % the 5th element is one
o_4 = [zeros(1,3) 1 zeros(1,2)]'; % the 6th element is one
z_3 = [ones(1,2) 0 ones(1,3)]';   % the 5th element is zero
z_4 = [ones(1,3) 0 ones(1,2)]';   % the 6th element is zero

for n = 1 : mod_sym_num/2
    dhat(2*n-1) = (HH1 * o_3)' * pinv(HH1 * diag(1 - (dd(2*n-1:2*n+4).^2).*z_3) * HH1' + eye(2*M,2*M)*Sigma) * (RR1((n-1)*M+1 : (n+1)*M) - HH1 * (dd(2*n-1:2*n+4) .* z_3)) + ...
                  (HH2 * o_3)' * pinv(HH2 * diag(1 - (dd(2*n-1:2*n+4).^2).*z_3) * HH2' + eye(2*M,2*M)*Sigma) * (RR2((n-1)*M+1 : (n+1)*M) - HH2 * (dd(2*n-1:2*n+4) .* z_3));
    dhat(2*n)   = (HH1 * o_4)' * pinv(HH1 * diag(1 - (dd(2*n-1:2*n+4).^2).*z_4) * HH1' + eye(2*M,2*M)*Sigma) * (RR1((n-1)*M+1 : (n+1)*M) - HH1 * (dd(2*n-1:2*n+4) .* z_4)) + ...
                  (HH2 * o_4)' * pinv(HH2 * diag(1 - (dd(2*n-1:2*n+4).^2).*z_4) * HH2' + eye(2*M,2*M)*Sigma) * (RR2((n-1)*M+1 : (n+1)*M) - HH2 * (dd(2*n-1:2*n+4) .* z_4));
%     dhat(2*n)   = (HH * o_6)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_6) * HH' + eye(3*M,3*M)*Sigma) * (RR((n-1)*M+1 : (n+2)*M) - HH * (dd(2*n-1:2*n+8) .* z_6));
    Mean(2*n-1) = (HH1 * o_3)' * pinv(HH1 * diag(1 - (dd(2*n-1:2*n+4).^2).*z_3) * HH1' + eye(2*M,2*M)*Sigma) * (HH1 * o_3) + ...
                  (HH2 * o_3)' * pinv(HH2 * diag(1 - (dd(2*n-1:2*n+4).^2).*z_3) * HH2' + eye(2*M,2*M)*Sigma) * (HH2 * o_3);
    Mean(2*n)   = (HH1 * o_4)' * pinv(HH1 * diag(1 - (dd(2*n-1:2*n+4).^2).*z_4) * HH1' + eye(2*M,2*M)*Sigma) * (HH1 * o_4) + ...
                  (HH2 * o_4)' * pinv(HH2 * diag(1 - (dd(2*n-1:2*n+4).^2).*z_4) * HH2' + eye(2*M,2*M)*Sigma) * (HH2 * o_4);
%     Mean(2*n)   = (HH * o_6)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_6) * HH' + eye(3*M,3*M)*Sigma) * (HH * o_6);
end

L_all = -2*real(dhat/2)./(1-real(Mean/2));

% for n = 1 : mod_sym_num/2
%     dhat(2*n-1) = (HH * o_5)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_5) * HH' + eye(3*M,3*M)*Sigma) * (RR((n-1)*M+1 : (n+2)*M) - HH * (dd(2*n-1:2*n+8) .* z_5));
%     dhat(2*n)   = (HH * o_6)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_6) * HH' + eye(3*M,3*M)*Sigma) * (RR((n-1)*M+1 : (n+2)*M) - HH * (dd(2*n-1:2*n+8) .* z_6));
%     Mean(2*n-1) = (HH * o_5)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_5) * HH' + eye(3*M,3*M)*Sigma) * (HH * o_5);
%     Mean(2*n)   = (HH * o_6)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_6) * HH' + eye(3*M,3*M)*Sigma) * (HH * o_6);
% end
% o_5 = [zeros(1,4) 1 zeros(1,5)]'; % the 5th element is one
% o_6 = [zeros(1,5) 1 zeros(1,4)]'; % the 6th element is one
% z_5 = [ones(1,4) 0 ones(1,5)]';   % the 5th element is zero
% z_6 = [ones(1,5) 0 ones(1,4)]';   % the 6th element is zero
% o_56 = [zeros(1,4) 1 1 zeros(1,4)]'; % the 5th 6th element is one
% z_56 = [ones(1,4) 0 0 ones(1,4)]';   % the 5th 6th element is zero

% for n = 1 : mod_sym_num/2
%     dhat(2*n-1) = (HH * o_5)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_56) * HH' + eye(3*M,3*M)*Sigma) * (RR((n-1)*M+1 : (n+2)*M) - HH * (dd(2*n-1:2*n+8) .* z_56));
%     dhat(2*n)   = (HH * o_6)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_56) * HH' + eye(3*M,3*M)*Sigma) * (RR((n-1)*M+1 : (n+2)*M) - HH * (dd(2*n-1:2*n+8) .* z_56));
%     Mean(2*n-1) = (HH * o_5)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_56) * HH' + eye(3*M,3*M)*Sigma) * (HH * o_5);
%     Mean(2*n)   = (HH * o_6)' * pinv(HH * diag(1 - (dd(2*n-1:2*n+8).^2).*z_56) * HH' + eye(3*M,3*M)*Sigma) * (HH * o_6);
% end

% L_all = -2*real(dhat)./(1-real(Mean));