diff --git a/.gitignore b/.gitignore
index 39108a0..ecda9b2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
 *.m~
+.DS_Store
 
diff --git a/jac.m b/jac.m
index 938c738..20a9610 100644
--- a/jac.m
+++ b/jac.m
@@ -18,47 +18,10 @@
 
 function C = jac(A,B)
 
-    mode = 0; % Exact Jacobian logarithm
-%	mode = 1; % Jacobian logarithm using a lookup table 
-%   mode = 2; % Approximate Jacobian logarithm
+    C = max(A,B) + log(1+exp(-abs(A-B)));
+%    C = max(A,B);
+    
+
 
-	if(A == -inf && B == -inf)
-        C = -inf;
-    else
-        if mode == 0
-        
-            C = max(A,B) + log(1+exp(-abs(A-B)));
-        
-        elseif mode == 1
-        
-            difference = abs(A-B);
-            if     difference >= 4.5
-                C = max(A,B);
-            elseif difference >= 2.252
-                C = max(A,B) + 0.05;
-            elseif difference >= 1.508
-                C = max(A,B) + 0.15;
-            elseif difference >= 1.05
-                C = max(A,B) + 0.25;
-            elseif difference >= 0.71
-                C = max(A,B) + 0.35;
-            elseif difference >= 0.433
-                C = max(A,B) + 0.45;
-            elseif difference >= 0.196
-                C = max(A,B) + 0.55;
-            else % difference >= 0
-                C = max(A,B) + 0.65;
-            end
-        
-        elseif mode == 2
-        
-            C = max(A,B);
-            
-        else
-            
-            error('Invalid Jacobian mode');
-            
-        end
-	end
 end
 	
\ No newline at end of file
diff --git a/main_inner.m b/main_inner.m
index 7d36ac3..9d8085d 100644
--- a/main_inner.m
+++ b/main_inner.m
@@ -22,7 +22,7 @@
 IA_count=11;
 
 % Channel SNR in dB
-SNR = -6;
+SNR = -4;
 
 % Noise variance
 N0 = 1/10^(SNR/10);
diff --git a/main_mod.m b/main_mod.m
new file mode 100644
index 0000000..eea9553
--- /dev/null
+++ b/main_mod.m
@@ -0,0 +1,88 @@
+% EXIT function for a soft-input soft-output demodulator
+% Copyright (C) 2008  Robert G. Maunder
+
+% This program is free software: you can redistribute it and/or modify it 
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your 
+% option) any later version.
+
+% This program is distributed in the hope that it will be useful, but 
+% WITHOUT ANY WARRANTY; without even the implied warranty of 
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
+% Public License for more details.
+
+% The GNU General Public License can be seen at http://www.gnu.org/licenses/.
+
+
+
+% Number of bits to encode
+bit_count=100000;
+
+% Number of a priori mutual informations to consider
+IA_count=11;
+
+% Channel SNR in dB
+SNR = 0;
+
+% Noise variance
+N0 = 1/10^(SNR/10);
+
+% A priori mutual informations to consider
+IA = 0.999*(0:1/(IA_count-1):1);
+
+% Initialise results
+IE_av=zeros(1,IA_count);
+IE_hist=zeros(1,IA_count);
+area=0.0;
+
+% Consider each a priori mutual information
+for IA_index = 1:IA_count
+
+    % Generate some random bits
+    bits  = round(rand(1,bit_count));
+
+    % Encode using a half-rate systematic recursive convolutional code having a single memory element
+    tx = modulate(bits);
+
+    % Rayleigh fading 
+    h = sqrt(1/2)*(randn(size(tx))+1i*randn(size(tx)));
+
+    % Noise
+    n = sqrt(N0/2)*(randn(size(tx))+1i*randn(size(tx)));
+
+    % Uncorrelated narrowband Rayleigh fading channel
+    rx = h.*tx + n;
+
+
+    % Generate the a priori LLRs having the a priori mutual information considered
+    apriori_llrs = generate_llrs(bits, IA(IA_index));
+  
+    % Do the BCJR
+    extrinsic_llrs = soft_demodulate(rx, h, N0, apriori_llrs);
+
+    % Measure the mutual information of the extrinsic LLRs
+    IE_hist(IA_index) = measure_mutual_information_histogram(extrinsic_llrs, bits);
+    IE_av(IA_index) = measure_mutual_information_averaging(extrinsic_llrs);
+
+    % Update the area beneath the EXIT function
+    if(IA_index > 1)
+       area = area + (IE_av(IA_index)+IE_av(IA_index-1))*(IA(IA_index)-IA(IA_index-1))/2;
+    end
+end
+
+
+% Plot EXIT function
+figure
+xlim([0 1]);
+ylim([0 1]);
+xlabel('Quality of input LLRs (a priori mutual information I_A)');
+ylabel('Quality of output LLRs (extrinsic mutual information I_E)');
+title(['EXIT function for SNR = ', num2str(SNR), ' dB']);
+hold on
+plot(IA,IE_hist,'r');
+plot(IA,IE_av,'b');
+legend({'True quality','Claimed quality'},'Location','northwest');
+
+% Display the area beneath the EXIT function
+annotation('textbox','String',{['Area = ', num2str(area)]},'LineStyle','none','Position',[0.7 0.1 0.2 0.1]);
+
diff --git a/modulate.m b/modulate.m
new file mode 100644
index 0000000..ef85529
--- /dev/null
+++ b/modulate.m
@@ -0,0 +1,40 @@
+% QPSK modulator using natural mapping
+% Copyright (C) 2010  Robert G. Maunder
+
+% This program is free software: you can redistribute it and/or modify it 
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your 
+% option) any later version.
+
+% This program is distributed in the hope that it will be useful, but 
+% WITHOUT ANY WARRANTY; without even the implied warranty of 
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
+% Public License for more details.
+
+% The GNU General Public License can be seen at http://www.gnu.org/licenses/.
+
+% bits is a k*n vector of bits
+% tx is a vector of complex symbols
+function tx = modulate(bits)
+
+    % Specify the constellation points and the bit mapping here
+    constellation_points = [+1+1i; -1+1i; -1-1i; +1-1i]/sqrt(2);
+    bit_labels = [0,0; 0,1; 1,1; 1,0];
+    
+    % Determine the number of bits per symbol and the number of constellation points here
+    k = size(bit_labels,2);
+    M = 2^k;
+    N = length(bits)/k;
+
+    
+    % Check that all the vectors and matrices have the correct dimensions
+    if ~isequal(size(constellation_points),[M,1]) || ~isequal(size(bit_labels),[M,k])
+        error('wrong dimensions');
+    end
+
+    symbols = bin2dec(num2str(reshape(bits,[k,N])'))'+1;
+    tx = constellation_points(symbols);
+    
+    
+end
+    
diff --git a/soft_demodulate.m b/soft_demodulate.m
new file mode 100644
index 0000000..462cf52
--- /dev/null
+++ b/soft_demodulate.m
@@ -0,0 +1,90 @@
+% Soft QPSK demodulator using natural mapping
+% Copyright (C) 2010  Robert G. Maunder
+
+% This program is free software: you can redistribute it and/or modify it 
+% under the terms of the GNU General Public License as published by the
+% Free Software Foundation, either version 3 of the License, or (at your 
+% option) any later version.
+
+% This program is distributed in the hope that it will be useful, but 
+% WITHOUT ANY WARRANTY; without even the implied warranty of 
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
+% Public License for more details.
+
+% The GNU General Public License can be seen at http://www.gnu.org/licenses/.
+
+% apriori_llrs is a 1xk vector of a priori LLRs
+% rx is a complex symbol
+% channel is a complex channel coefficient
+% N0 is the noise power spectral density
+% extrinsic_llrs is a 1xk vector of extrinsic LLRs
+function extrinsic_llrs = soft_demodulate(rx, channel, N0, apriori_llrs)
+
+    % Specify the constellation points and the bit mapping here
+    constellation_points = [+1+1i; -1+1i; -1-1i; +1-1i]/sqrt(2);
+    bit_labels = [0,0; 0,1; 1,0; 1,1];
+    
+    % Determine the number of bits per symbol and the number of constellation points here
+    k = size(bit_labels,2);
+    M = 2^k;
+    N = length(rx);
+
+    % Check that all the vectors and matrices have the correct dimensions
+    if ~isequal(size(constellation_points),[M,1]) || ~isequal(size(bit_labels),[M,k])
+        error('wrong dimensions');
+    end
+    
+    if length(channel) ~= length(rx) && length(channel) ~= 1
+        error('wrong dimensions');
+    end
+    
+    
+    
+aposteriori_symbol_LLRs = zeros(M,N);
+
+% Put the influence of the received signals into the symbol LLRs
+for perm_index = 1:M
+    aposteriori_symbol_LLRs(perm_index,:) = -abs(rx-channel*constellation_points(perm_index)).^2./N0;
+end
+
+if exist('apriori_llrs','var')
+    % Put the influence of the apriori LLRs into the symbol LLRs
+    for bit_index = 1:k       
+%        aposteriori_symbol_LLRs(:,bit_permutations(bit_index,:) == 0) = aposteriori_symbol_LLRs(:,bit_permutations(bit_index,:) == 0) + repmat(apriori_llrs(bit_index:bits_per_symbol:end),[1,permutations/2]);
+        
+        for perm_index = 1:M
+            if bit_labels(perm_index,bit_index) == 0
+                aposteriori_symbol_LLRs(perm_index,:) = aposteriori_symbol_LLRs(perm_index,:) + apriori_llrs(bit_index:k:end);
+            end
+        end   
+    end
+end
+
+% Extract the aposteriori LLRs from the symbol LLRs
+aposteriori_llrs = zeros(1,N*k);
+for bit_index = 1:k
+    p0 = -inf(1,N);
+    p1 = -inf(1,N);
+
+    for perm_index = 1:M
+        if bit_labels(perm_index,bit_index) == 0
+            p0 = jac(p0, aposteriori_symbol_LLRs(perm_index,:));
+        else
+            p1 = jac(p1, aposteriori_symbol_LLRs(perm_index,:));
+        end
+    end   
+    
+    aposteriori_llrs(bit_index:k:end) = p0-p1;
+end
+        
+if exist('apriori_llrs','var')
+    % Remove the apriori from the aposteriori to get the extrinsic
+    extrinsic_llrs = aposteriori_llrs - apriori_llrs;
+else
+    extrinsic_llrs = aposteriori_llrs;
+end
+
+
+    
+ 
+end