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);
d_r = zeros(M, mod_sym_num);
dhat = zeros(1, mod_sym_num);
Mean = zeros(1, mod_sym_num);
j = sqrt(-1);

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
        d_r(l,2*k-1) = real(conj(H(l,1))*R(l,2*k-1)/j + H(l,2)*conj(R(l,2*k)));     % d(2n),  n=0,1,2...
        d_r(l,2*k) = real(H(l,1)*conj(R(l,2*k)) - conj(H(l,2))*R(l,2*k-1)/j);       % d(2n+1),n=0,1,2...
    end
end

d_r = [zeros(M,2) d_r zeros(M,2)];
RR = reshape(d_r,M*(mod_sym_num+4),1);

L_a = [zeros(1,4) L_a zeros(1,4)]'; % 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 = sum(abs(H(1,1))^2)+sigma^2;

for l = 1 : M
    h(l,:) = real([-j*conj(H(l,1))*H(l,2)*c(2); 0; ((abs(H(l,1))^2+abs(H(l,2))^2)*c(1)-2*j*conj(H(l,1))*H(l,2)*c(2)); 0; -j*conj(H(l,1))*H(l,2)*c(2)])/sqrt(2*M); % Channel Coefficients
end

HH = [h zeros(M,4); zeros(M,1) h zeros(M,3); zeros(M,2) h zeros(M,2); zeros(M,3) h zeros(M,1); zeros(M,4) h]; % Coefficients Matrix

o_5 = [zeros(1,4) 1 zeros(1,4)]'; % the 5th element is one
z_5 = [ ones(1,4) 0  ones(1,4)]'; % the 5th element is zero

for n = 1 : mod_sym_num
    dhat(n) = (HH * o_5)' * pinv(HH * diag(1 - (dd(n:n+8).^2).*z_5) * HH' + eye(M*5,M*5)*Sigma) * (RR((n-1)*M+1 : (n+4)*M) - HH * (dd(n:n+8).* z_5));
    Mean(n) = (HH * o_5)' * pinv(HH * diag(1 - (dd(n:n+8).^2).*z_5) * HH' + eye(M*5,M*5)*Sigma) * (HH * o_5);
end

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