function L_all = mmse( 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);
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
% R = [zeros(M,1) R zeros(M,1)];
% RR = reshape(R,M*(mod_sym_num+2),1);
% rr = [zeros(M,1) R zeros(M,1)];
% RR = reshape(rr,M*(mod_sym_num+2),1);
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);

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 = 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); % 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); % ISI coefficients
end
% HH = [h2 zeros(M,2);h1 zeros(M,2);zeros(M,2) h2;zeros(M,2) h1];
HH1 = [h2 zeros(M,2);h1 zeros(M,2)];
HH2 = [zeros(M,2) h2;zeros(M,2) h1];
HH =  [HH1;HH2];

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) = (HH * o_3)'*(pinv(HH * diag(1 - (dd(2*n-1:2*n+4).^2).*z_3) * HH' + eye(4*M,4*M)*Sigma)) * ( [RR2((n-1)*M+1 : (n)*M);RR1((n-1)*M+1 : (n)*M);RR2((n)*M+1 : (n+1)*M);RR1((n)*M+1 : (n+1)*M)] - HH * (dd(2*n-1:2*n+4) .* z_3));
     dhat(2*n)   = (HH * o_4)'*(pinv(HH * diag(1 - (dd(2*n-1:2*n+4).^2).*z_4) * HH' + eye(4*M,4*M)*Sigma)) * ( [RR2((n-1)*M+1 : (n)*M);RR1((n-1)*M+1 : (n)*M);RR2((n)*M+1 : (n+1)*M);RR1((n)*M+1 : (n+1)*M)] - HH * (dd(2*n-1:2*n+4) .* z_4));
     Mean(2*n-1) = (HH * o_3)'*(pinv(HH * diag(1 - (dd(2*n-1:2*n+4).^2).*z_3) * HH' + eye(4*M,4*M)*Sigma)) * (HH * o_3);
     Mean(2*n)   = (HH * o_4)'*(pinv(HH * diag(1 - (dd(2*n-1:2*n+4).^2).*z_4) * HH' + eye(4*M,4*M)*Sigma)) * (HH * o_4);
end

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