msksim.f90 works, but needs more work.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6493 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

lib/ldpc/Makefile.msk_osx

@@ -27,7 +27,7 @@ all:	msksim
 
 OBJS = msksim.o alloc.o rcode.o channel.o dec.o enc.o \
 	intio.o blockio.o check.o open.o mod2dense.o \
-	mod2sparse.o mod2convert.o distrib.o rand.o
+	mod2sparse.o mod2convert.o distrib.o rand.o gran.o
 
 msksim:	$(OBJS)
 	$(FC) -o msksim $(OBJS) 

lib/ldpc/dec.c

@@ -384,19 +384,23 @@ void iterprp
   }
 }
 
-void ldpc_decode_ ( double *lratio, char *dblk, char *pchk, double *bprb )
+void ldpc_decode_ ( double lratio[], char decoded[], int *max_iterations, int *niterations)
 {
-  int i, n;
-  n = prprp_decode ( H, lratio, dblk, pchk, bprb );
-  printf("in ldpc_decode n=%d\n");
-  printf("dblk: ");
-  for( i=0; i<128; i++) {
-    printf("%d",dblk[i]);
+  int i, j, valid;
+  char dblk[N],pchk[M];
+  double bprb[N];
+    
+  max_iter=*max_iterations;
+  *niterations = prprp_decode ( H, lratio, dblk, pchk, bprb );
+  valid = check( H, dblk, pchk )==0;
+  if( !valid ) {
+      *niterations=-1;
+  };
+
+  j=0;
+  for( i=M; i<N; i++ ) {
+    decoded[j]=dblk[cols[i]];
+    j=j+1;
   }
-  printf("\n");
-  printf("pchk: ");
-  for( i=0; i<46; i++) {
-    printf("%d",pchk[i]);
-  }
-  printf("\n");
+
 } 

lib/ldpc/enc.c

@@ -159,11 +159,20 @@ char msg[82],
 char cdw[128]
 ){
   int i;
+  char checks[46];
+
   printf("msg: ");
   for (i=0; i<82; i++) {
   printf("%d", msg[i]); }
   printf("\n");
   sparse_encode(msg,cdw);
+  mod2sparse_mulvec (H, cdw, checks);
+  for (i = 0; i<46; i++) {
+    if( checks[i] == 1 ) {
+      printf("Failed to encode.\n");
+      abort();
+    }
+  } 
   printf("cdw: ");
   for (i=0; i<128; i++) {
   printf("%d", cdw[i]); }

lib/ldpc/msksim.f90

@@ -5,25 +5,60 @@ use, intrinsic :: iso_c_binding
 parameter (N=128, M=46, K=82) ! M and N are global variables on the C side. 
 integer(1) message(1:K)
 integer(1) codeword(1:N)
-integer(1) decoded(1:N)
-integer(1) pchk(1:M)
-real*8 lratio(N), bitprobs(N)
+integer(1) decoded(1:K)
+real*8 lratio(N)
+character(50) pchk_file,gen_file
 
-write(*,*) "calling init_ldpc"
-call init_ldpc()
+pchk_file="./jtmode_codes/ldpc-128-82-sf11.pchk"
+gen_file="./jtmode_codes/ldpc-128-82-sf11.gen"
 
-message(1:K)=0
-message(10)=1
-write(*,*) "calling ldpc_encode"
+call init_ldpc(trim(pchk_file)//char(0),trim(gen_file)//char(0))
+
+
+message(1:41)=1
+message(42:82)=0
 call ldpc_encode(message,codeword)
-write(*,*) "calling ldpc_decode"
-do i=1,N
-lratio(i)=exp(2.0*(codeword(i)-0.5))
-enddo
-lratio(10)=10.0
-bitprobs(1:N)=0.0 ! shouldn't need this
-call ldpc_decode(lratio, decoded, pchk, bitprobs)
-do i=1,N
-write(*,*) i,bitprobs(i), lratio(i)
+
+max_iterations=10
+ntrials=1000000
+rate=82.0/128.0
+
+write(*,*) "Eb/N0   ngood    nundetected"
+do idb = 0, 11
+  db=idb/2.0-0.5
+  sigma=1/sqrt( 2*rate*(10**(db/10.0)) )
+
+  ngood=0
+  nue=0
+  aviter=0.0
+
+  do itrial=1, ntrials
+
+    do i=1,N
+      rxdata = 2.0*(codeword(i)-0.5) + sigma*gran()
+      lratio(i)=exp(2.0*rxdata/(sigma*sigma))
+    enddo
+
+    call ldpc_decode(lratio, decoded, max_iterations, niterations)
+
+    if( niterations .ge. 0 ) then
+      nueflag=0
+      do i=1,K
+        if( message(i) .ne. decoded(i) ) then
+          nueflag=1
+        endif
+      enddo
+      if( nueflag .eq. 1 ) then
+        nue=nue+1
+      else
+        ngood=ngood+1;
+      endif
+    endif
+
+  enddo
+
+  write(*,"(f4.1,1x,i8,1x,i8)") db,ngood,nue
+
 enddo
+
 end program msksim

lib/ldpc/rcode.c

@@ -199,19 +199,11 @@ garbled:
   exit(1);
 }
 
-// Fortran interface routine
-void init_ldpc_ ( )
+// Fortran interface routine for WSJT-X
+void init_ldpc_ (char *pfile, char *gfile )
 {
     char *pchk_file,*gen_file;
-    pchk_file="./jtmode_codes/ldpc-128-82-sf11.pchk";
-    gen_file="./jtmode_codes/ldpc-128-82-sf11.gen";
- 
-    printf("pchk_file %s\n",pchk_file);
-    printf("gen_file %s\n",gen_file);
-
-    read_pchk( pchk_file );
-    read_gen( gen_file, 0, 0 );
-
-    printf("N %d M %d\n",N,M);
+    read_pchk( pfile );
+    read_gen( gfile, 0, 0 );
 }