Comments welcomed to progr. style / speed of my first ocaml program "Error correcting LDPC decoder"

[Andries Hekstra <andries.hekstra@philips.com>]

[-- Attachment #1.1: Type: text/plain, Size: 1277 bytes --]

Dear All,

I am an applied scientist in the areas of signal processing and error 
correction at Philips Research. Having done the programming work of my 
computer simulations in C/C++ for almost 10-15 years with increasingly 
mixed feelings about the languange, a friend pointed me to functional 
programming. This way I found ocaml and I must say I am very enthusiastic 
about it. The built-in arrays, the functions that are really meta 
functions and remove a lot of tedious programming and debugging. Also the 
combined availability of an interpreter and compiler is comforbable and 
not unlike matlab (which I hate). 

I am also intrigued by the promise of functional languages when multi-core 
computers will be the defacto standard. 

Attached please find the first program I wrote. It has about the same 
speed as its C++ counterpart (slightly faster). I welcome comments w.r.t. 
programming style, esp. when it affects the speed of the program.

Regards,

Andries Hekstra




------------------------------------------------------------------------
Dr. Ir. Andries P. Hekstra
Philips Research 
High Tech Campus 37 
5656 AG Eindhoven
Tel./Fax/Secr. +31 40 27 42048/42566/44051 

If you need anti-RSI break software : Look at http://www.workrave.org 
(good and for free)

[-- Attachment #1.2: Type: text/html, Size: 1628 bytes --]

[-- Attachment #2: Main.ml --]
[-- Type: application/octet-stream, Size: 16249 bytes --]

(* Parameters of the LDPC code, left and right degree*)
let    ldpc_j = 3
let    ldpc_k = 30
let    codeRate = 1. -. float_of_int ldpc_j /. float_of_int ldpc_k
let    targetCodewordLength = (ldpc_k * 300)
let    blockSize = targetCodewordLength / ldpc_k
let    codewordLength = blockSize * ldpc_k
let    numberOfParityChecks = blockSize * ldpc_j
let    fudgeFactor = 0.75
let    maxNumberOfIterations = 30
;;

print_string "LDPC decoder for rate 1 - ";
print_int ldpc_j;
print_string "/";
print_int ldpc_k;
print_string " Gallager code with length ";
print_int codewordLength;
flush stdout
;;

(* Construction of a sparse parity check matrix*)
let height a = Array.length a
let width a = Array.length a.(0)

let applyToMatrix f a =
    let rec applyToMatrixAux i j f a = 
        if (i >= 0) then begin 
            a.(i).(j) <- f i j a.(i).(j); 
            if (j-1 >= 0) 
                then applyToMatrixAux i     (j-1)         f a 
                else applyToMatrixAux (i-1) (width a - 1) f a
        end
        in applyToMatrixAux (height a - 1) (width a - 1) f a

let order =       Array.create_matrix ldpc_j codewordLength 0
let permutation = Array.create_matrix ldpc_j codewordLength 0
;;

applyToMatrix (fun j -> fun i -> fun x -> (Random.int (100 * codewordLength))) order;; 
applyToMatrix (fun j -> fun i -> fun x -> i) permutation;;

let swapElement i j a = let x = a.(i) in a.(i) <- a.(j); a.(j) <- x

let qsort op a = 
    let l = Array.length a in     
    let rec qsortAux i j op a =
        if (i < j) then begin 
            let x = a.((i + j) / 2) in
            let i1 = ref i in
            let j1 = ref j in                                
            
            while (!i1 <= !j1) do 
                while (op a.(!i1) x) do incr i1 done;
                while (op x a.(!j1)) do decr j1 done;
                if (i1 <= j1) then begin
                    swapElement !i1 !j1 a; 
                    incr i1;
                    decr j1
                end
            done;                                        
            
            qsortAux !i1 j op a;
            qsortAux i !j1 op a;                                 
         end                
        in qsortAux 0 (l - 1) op a;;

let qsort2 op a b = 
    let l = Array.length a in     
    let rec qsort2Aux i j op a b =
        if (i < j) then begin 
            let x = a.((i + j) / 2) in
            let i1 = ref i in
            let j1 = ref j in
                            
            while (!i1 <= !j1) do 
                while (op a.(!i1) x) do incr i1 done;
                while (op x a.(!j1)) do decr j1 done;
                if (i1 <= j1) then begin
                    swapElement !i1 !j1 a; 
                    swapElement !i1 !j1 b; 
                    incr i1;
                    decr j1
                 end
            done;                            
            qsort2Aux !i1 j op a b;
            qsort2Aux i !j1 op a b;        
        end                
        in qsort2Aux 0 (l - 1) op a b;;

let sortRows a = 
    let rec sortRowsAux i a = if (i >= 0) then begin qsort (<) a.(i); sortRowsAux (i-1) a end  
        in sortRowsAux (Array.length a - 1) a

let sortRows2 a b = 
    let rec sortRows2Aux i a b = if (i >= 0) then begin qsort2 (<) a.(i) b.(i); sortRows2Aux (i-1) a b end  
        in sortRows2Aux (Array.length a - 1) a b;;

sortRows2 order permutation;;                                                                                

let inversePermutation = Array.create_matrix ldpc_j codewordLength 0;;

for j = 0 to ldpc_j - 1 do
    for i = 0 to codewordLength - 1 do 
        inversePermutation.(j).(permutation.(j).(i)) <- i
    done
done;;     

let oneColumnInRow = Array.create_matrix numberOfParityChecks ldpc_k 0;;
let jOneColumnInRow = Array.create_matrix numberOfParityChecks ldpc_k 0;;

    (* Here, i indexes parity checks *)
applyToMatrix (fun i -> fun k -> fun x -> inversePermutation.(i / blockSize).(k + ldpc_k * (i mod blockSize))) oneColumnInRow;;
applyToMatrix (fun i -> fun k -> fun x -> i / blockSize) jOneColumnInRow;;

sortRows2 oneColumnInRow jOneColumnInRow;;

let oneRowInColumn = Array.create_matrix codewordLength ldpc_j 0;;
let kOneRowInColumn = Array.create_matrix codewordLength ldpc_j 0;;

for j = 0 to ldpc_j-1 do
    for r = 0 to blockSize -1 do 
        let i = r + j * blockSize in 
        for k = 0 to ldpc_k-1 do
            oneRowInColumn .(oneColumnInRow.(i).(k)).(j) <- i;
            kOneRowInColumn.(oneColumnInRow.(i).(k)).(j) <- k;
        done                        
    done
done;;

sortRows2 oneRowInColumn kOneRowInColumn;;

let numberOfColumnPairsWithOverlap = ref 0;; 
for i1 = 0 to codewordLength - 1 do
    for i2 = 0 to i1-1 do
        let overlap = ref 0 in begin 
            for j1 = 0 to ldpc_j-1 do
                for j2 = 0 to ldpc_j-1 do
                    if (oneRowInColumn.(i1).(j1) = oneRowInColumn.(i2).(j2)) then incr overlap;
                done
            done;
            if (!overlap > 1) then incr numberOfColumnPairsWithOverlap;
        end        
    done
done;;
print_string "\nnumberOfColumnPairsWithOverlap = ";;
print_int !numberOfColumnPairsWithOverlap;;

(* Random bit vector *)
let bitArrayFromInt n iii = 
    let rec bitArrayFromIntAux i n iii = 
        if (i < n) then 
            Array.append [| (iii mod 2) |]  (bitArrayFromIntAux (i+1) n (iii / 2)) 
        else [||]
        in bitArrayFromIntAux 0 n iii 

let makeRandomBit a = 
    let rec makeRandomBitAux i a = 
        if (i > 0) then begin
            let b = bitArrayFromInt 8 (Random.int 256) in
            let i1 = ref i in 
            let count = ref 0 in begin       
                while (!i1 >= 0) & (!count < 8) do
                    a.(!i1) <- b.(!count);
                    decr i1;
                    incr count;
                done;
                makeRandomBitAux !i1 a
            end    
        end         
    in makeRandomBitAux (Array.length a - 1) a

(* Normal distribution *)    
let randomGaussian sigma =
    let u1 = Random.float 1. in
    let u2 = Random.float 1. in        
    let r = sqrt (-2. *. log (u1)) in
    let pi = 2. *. asin 1. in
    let phi = 2. *. pi *. u2 
        in sigma *. r *. cos (phi)        
         
(* Message passing functions *)
let sum a = let rec sum i a = if (i > 0) then a.(i) +. sum (i-1) a else if (i = 0) then a.(i) else 0.
    in sum ((Array.length a) - 1) a
let hardBit softBit = if (softBit >= 0.) then 0 else 1
let reliability softBit = if (softBit >= 0.) then softBit else (-. softBit)

(* LDPC decode *)
let softBitOut channelLLR1 messageSymbolNodeIn1 = channelLLR1 +. sum messageSymbolNodeIn1 
let hardBitOut channelLLR1 messageSymbolNodeIn1 = hardBit (softBitOut channelLLR1 messageSymbolNodeIn1)         

let processCheckNode channelLLR 
                     oneColumnInRow 
                     jOneColumnInRow
                     syndrome  
                     messagesSymbolNodeIn 
                     messageCheckNodeIn 
                     messageCheckNodeOut = 
                     
    let messageSymbolNodeOut channelLLR1 messageSymbolNodeIn1 j = 
            (softBitOut channelLLR1 messageSymbolNodeIn1) -. messageSymbolNodeIn1.(j) in begin
            
        (* Get the input messages to the current check node *)
        for k = 0 to ldpc_k-1 do 
            messageCheckNodeIn.(k) <- fudgeFactor *. 
                                          (messageSymbolNodeOut channelLLR.(oneColumnInRow.(k)) 
                                                                messagesSymbolNodeIn.(oneColumnInRow.(k)) 
                                                                jOneColumnInRow.(k))                                    
        done; 

        (* Compute the total EXOR, the min. and second min. reliability and the index *)
        let k1 = ref 0 in  
        let minReliability1 = ref 1e10 in 
        let minReliability2 = ref 1e10 in
        let totalExor = ref syndrome in begin                
           
            for k = 0 to ldpc_k-1 do 
                totalExor := !totalExor + hardBit messageCheckNodeIn.(k);
                
                let h = reliability messageCheckNodeIn.(k) in  
                    if !minReliability1 >= h then begin
                        minReliability2 := !minReliability1;
                        minReliability1 := h;
                        k1 := k;
                    end else begin
                        if !minReliability2 >= h then minReliability2 := h                                
                    end
            done; 
        
            totalExor := !totalExor mod 2;          
            
            (* Compute the output messages *)
            for k = 0 to ldpc_k-1 do 
                let h = if k = !k1 then !minReliability2 else !minReliability1 in
                    if ((!totalExor + (hardBit messageCheckNodeIn.(k))) mod 2) = 0 then
                        messageCheckNodeOut.(k) <- h
                    else 
                        messageCheckNodeOut.(k) <- (-. h)
            done;
            
            (* Assign the output messages *)
            for k = 0 to ldpc_k-1 do 
                messagesSymbolNodeIn.(oneColumnInRow.(k)).(jOneColumnInRow.(k)) <- messageCheckNodeOut.(k); 
            done;
            
            (!totalExor = 0)    
        end
    end 
          
let ldpcDecode transmittedCodeword 
               syndrome
               channelLLR 
               decodedCodeword 
               oneColumnInRow 
               jOneColumnInRow 
               messagesSymbolNodeIn 
               numberOfIterations =
    
    let messageCheckNodeIn = Array.create ldpc_k 0. in
    let messageCheckNodeOut = Array.create ldpc_k 0. in
    let allZeroesSyndrome = ref false in begin
        numberOfIterations := 0;
    
        while (!numberOfIterations < maxNumberOfIterations) & not !allZeroesSyndrome do 
            allZeroesSyndrome := true;
            for row = 0 to numberOfParityChecks - 1 do 
                if not (processCheckNode channelLLR 
                                         oneColumnInRow.(row) 
                                         jOneColumnInRow.(row) 
                                         syndrome.(row)
                                         messagesSymbolNodeIn
                                         messageCheckNodeIn
                                         messageCheckNodeOut) then
                    allZeroesSyndrome := false;                    
            done;    
            incr numberOfIterations;
        done;  
        
        let numberOfBitErrors = ref 0 in begin
            for column = 0 to codewordLength - 1 do
                decodedCodeword.(column) <- hardBitOut channelLLR.(column) messagesSymbolNodeIn.(column);
                if (decodedCodeword.(column) <> transmittedCodeword.(column)) then
                    incr numberOfBitErrors;    
            done;
    
            !numberOfBitErrors;   
        end
    end
        
(* Simulation parameters *)
let ebOverNoDbLow = 3.4
let ebOverNoDbHigh = 4.
(*let ebOverNoDbLow = 10.0
let ebOverNoDbHigh = 10.0*)
let ebOverNoDbStep = 0.1
let numberOfCodewordsToSimulate = 30000

let ebOverNoDb = ref ebOverNoDbLow
let totalNumberOfFrameErrors = ref 1  
;;

while (!ebOverNoDb <= ebOverNoDbHigh) & (!totalNumberOfFrameErrors > 0) do

    print_string "\n\nEbOverNoDb = ";
    print_float  !ebOverNoDb;    
    print_string ", time = ";
    print_float (Sys.time());
    print_string " seconds!\n";
    flush stdout;

    let ebOverNo = 10. ** (!ebOverNoDb /. 10.) in 
    let lc = 4. *. ebOverNo *. codeRate in 
    let sigmaNoise = 1. /. sqrt (2. *. ebOverNo *. codeRate) in 
    
    let totalNumberOfChannelBitErrors = ref 0. in
    let totalNumberOfBits = ref 0. in
    let totalNumberOfDecodedBitErrors = ref 0. in 
    let totalNumberOfFrames = ref 0 in 
    let totalNumberOfIterations = ref 0. in 
        
    let transmittedCodeword = Array.create codewordLength 0 in 
    let decodedCodeword = Array.create codewordLength 0 in 
    let syndrome = Array.create numberOfParityChecks 0 in 
    let channelLLR = Array.create codewordLength 0. in 
    
    let messagesSymbolNodeIn = Array.create_matrix codewordLength ldpc_j 0. in begin
    
        totalNumberOfFrameErrors := 0; 
    
        while (!totalNumberOfFrames < numberOfCodewordsToSimulate) 
                & ((!totalNumberOfDecodedBitErrors < 500.) or (!totalNumberOfFrameErrors < 200)) do    
        
            makeRandomBit transmittedCodeword;
                                
            for r = 0 to numberOfParityChecks-1 do
                syndrome.(r) <- 0;
                for k = 0 to ldpc_k-1 do
                    syndrome.(r) <- ((syndrome.(r) + transmittedCodeword.(oneColumnInRow.(r).(k))) mod 2)
                done
            done;
            
            for i = 0 to codewordLength - 1 do
                channelLLR.(i) <- float_of_int (1 - 2*transmittedCodeword.(i)) +. randomGaussian sigmaNoise
            done;
            
            let numberOfChannelBitErrors = ref 0 in begin
                for i = 0 to codewordLength-1 do
                    if transmittedCodeword.(i) <> hardBit channelLLR.(i) then incr numberOfChannelBitErrors
                done;                
                
                totalNumberOfChannelBitErrors := !totalNumberOfChannelBitErrors +. float_of_int !numberOfChannelBitErrors;
                totalNumberOfBits := !totalNumberOfBits +. float_of_int codewordLength;
            end;         
        
            let numberOfIterations = ref 0 in  
                let numberOfDecodedBitErrors = ldpcDecode transmittedCodeword 
                                                          syndrome
                                                          channelLLR 
                                                          decodedCodeword 
                                                          oneColumnInRow 
                                                          jOneColumnInRow 
                                                          messagesSymbolNodeIn 
                                                          numberOfIterations in begin
                    totalNumberOfDecodedBitErrors := !totalNumberOfDecodedBitErrors +. float_of_int numberOfDecodedBitErrors; 
                    totalNumberOfIterations := !totalNumberOfIterations +. float_of_int !numberOfIterations;
                    if numberOfDecodedBitErrors > 0 then incr totalNumberOfFrameErrors;
                    incr totalNumberOfFrames;
                end;

            if (!totalNumberOfFrames mod 100) = 0 then begin                       
                let avgChannelBER = !totalNumberOfChannelBitErrors /. !totalNumberOfBits in 
                let avgDecodedBER = !totalNumberOfDecodedBitErrors /. !totalNumberOfBits in 
                let avgFER = float_of_int !totalNumberOfFrameErrors /. float_of_int !totalNumberOfFrames in 
                let avgIter = !totalNumberOfIterations /. float_of_int !totalNumberOfFrames in 
                begin
                    print_int !totalNumberOfFrames;
                    print_string "\tChannelBER = ";
                    print_float avgChannelBER;
                    print_string "\tDecodedBER = ";
                    print_float avgDecodedBER;
                    print_string "\tFER = ";
                    print_float avgDecodedBER;
                    print_string "\tavgIter = ";
                    print_float avgIter;
                    print_string "\n";
                    flush stdout;
                end     
            end
        done
    end;
    ebOverNoDb := !ebOverNoDb +. ebOverNoDbStep;
done;;
            
print_string "\n\nSuccesful completion in ";;
print_float (Sys.time());;
print_string " seconds!\n";;