diff --git a/libm65/CMakeLists.txt b/libm65/CMakeLists.txt
index f80199df4..77cd2454e 100644
--- a/libm65/CMakeLists.txt
+++ b/libm65/CMakeLists.txt
@@ -71,6 +71,7 @@ set (FSRCS
   ftnquit.f90
   gen65.f90
   getdphi.f90
+  graycode65.f90
   iqcal.f90
   iqfix.f90
   jt65code.f90
@@ -138,6 +139,7 @@ set (FSRCS
 set (CSRCS
   decode_rs.c
   encode_rs.c
+  ftrsd2.c
   gran.c
   igray.c
   init_rs.c
diff --git a/libm65/decode_rs.c b/libm65/decode_rs.c
index aecc47faf..91f582ac1 100644
--- a/libm65/decode_rs.c
+++ b/libm65/decode_rs.c
@@ -1,6 +1,7 @@
 /* Reed-Solomon decoder
  * Copyright 2002 Phil Karn, KA9Q
  * May be used under the terms of the GNU General Public License (GPL)
+ * Modified by Steve Franke, K9AN, for use in a soft-symbol RS decoder
  */
 
 #ifdef DEBUG
@@ -21,243 +22,247 @@
 #endif
 
 int DECODE_RS(
-#ifdef FIXED
-DTYPE *data, int *eras_pos, int no_eras,int pad){
-#else
-void *p,DTYPE *data, int *eras_pos, int no_eras){
-  struct rs *rs = (struct rs *)p;
+#ifndef FIXED
+              void *p,
 #endif
-  int deg_lambda, el, deg_omega;
-  int i, j, r,k;
-  DTYPE u,q,tmp,num1,num2,den,discr_r;
-  DTYPE lambda[NROOTS+1], s[NROOTS];	/* Err+Eras Locator poly
-					 * and syndrome poly */
-  DTYPE b[NROOTS+1], t[NROOTS+1], omega[NROOTS+1];
-  DTYPE root[NROOTS], reg[NROOTS+1], loc[NROOTS];
-  int syn_error, count;
-
-#ifdef FIXED
-  /* Check pad parameter for validity */
-  if(pad < 0 || pad >= NN)
-    return -1;
-#endif
-
-  /* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
-  for(i=0;i<NROOTS;i++)
-    s[i] = data[0];
-
-  for(j=1;j<NN-PAD;j++){
-    for(i=0;i<NROOTS;i++){
-      if(s[i] == 0){
-	s[i] = data[j];
-      } else {
-	s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR+i)*PRIM)];
-      }
-    }
-  }
-
-  /* Convert syndromes to index form, checking for nonzero condition */
-  syn_error = 0;
-  for(i=0;i<NROOTS;i++){
-    syn_error |= s[i];
-    s[i] = INDEX_OF[s[i]];
-  }
-
-  if (!syn_error) {
-    /* if syndrome is zero, data[] is a codeword and there are no
-     * errors to correct. So return data[] unmodified
-     */
-    count = 0;
-    goto finish;
-  }
-  memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
-  lambda[0] = 1;
-
-  if (no_eras > 0) {
-    /* Init lambda to be the erasure locator polynomial */
-    lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))];
-    for (i = 1; i < no_eras; i++) {
-      u = MODNN(PRIM*(NN-1-eras_pos[i]));
-      for (j = i+1; j > 0; j--) {
-	tmp = INDEX_OF[lambda[j - 1]];
-	if(tmp != A0)
-	  lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
-      }
-    }
-
-#if DEBUG >= 1
-    /* Test code that verifies the erasure locator polynomial just constructed
-       Needed only for decoder debugging. */
+              DTYPE *data, int *eras_pos, int no_eras, int calc_syn){
     
-    /* find roots of the erasure location polynomial */
-    for(i=1;i<=no_eras;i++)
-      reg[i] = INDEX_OF[lambda[i]];
-
-    count = 0;
-    for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
-      q = 1;
-      for (j = 1; j <= no_eras; j++)
-	if (reg[j] != A0) {
-	  reg[j] = MODNN(reg[j] + j);
-	  q ^= ALPHA_TO[reg[j]];
-	}
-      if (q != 0)
-	continue;
-      /* store root and error location number indices */
-      root[count] = i;
-      loc[count] = k;
-      count++;
-    }
-    if (count != no_eras) {
-      printf("count = %d no_eras = %d\n lambda(x) is WRONG\n",count,no_eras);
-      count = -1;
-      goto finish;
+#ifndef FIXED
+    struct rs *rs = (struct rs *)p;
+#endif
+    int deg_lambda, el, deg_omega;
+    int i, j, r,k;
+    DTYPE u,q,tmp,num1,num2,den,discr_r;
+    DTYPE lambda[NROOTS+1];	// Err+Eras Locator poly
+    static DTYPE s[51];					 // and syndrome poly
+    DTYPE b[NROOTS+1], t[NROOTS+1], omega[NROOTS+1];
+    DTYPE root[NROOTS], reg[NROOTS+1], loc[NROOTS];
+    int syn_error, count;
+    
+    if( calc_syn ) {
+        /* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
+        for(i=0;i<NROOTS;i++)
+            s[i] = data[0];
+        
+        for(j=1;j<NN;j++){
+            for(i=0;i<NROOTS;i++){
+                if(s[i] == 0){
+                    s[i] = data[j];
+                } else {
+                    s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR+i)*PRIM)];
+                }
+            }
+        }
+        
+        /* Convert syndromes to index form, checking for nonzero condition */
+        syn_error = 0;
+        for(i=0;i<NROOTS;i++){
+            syn_error |= s[i];
+            s[i] = INDEX_OF[s[i]];
+        }
+        
+        
+        if (!syn_error) {
+            /* if syndrome is zero, data[] is a codeword and there are no
+             * errors to correct. So return data[] unmodified
+             */
+            count = 0;
+            goto finish;
+        }
+        
     }
+    
+    memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
+    lambda[0] = 1;
+    
+    if (no_eras > 0) {
+        /* Init lambda to be the erasure locator polynomial */
+        lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))];
+        for (i = 1; i < no_eras; i++) {
+            u = MODNN(PRIM*(NN-1-eras_pos[i]));
+            for (j = i+1; j > 0; j--) {
+                tmp = INDEX_OF[lambda[j - 1]];
+                if(tmp != A0)
+                    lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
+            }
+        }
+        
+#if DEBUG >= 1
+        /* Test code that verifies the erasure locator polynomial just constructed
+         Needed only for decoder debugging. */
+        
+        /* find roots of the erasure location polynomial */
+        for(i=1;i<=no_eras;i++)
+            reg[i] = INDEX_OF[lambda[i]];
+        
+        count = 0;
+        for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
+            q = 1;
+            for (j = 1; j <= no_eras; j++)
+                if (reg[j] != A0) {
+                    reg[j] = MODNN(reg[j] + j);
+                    q ^= ALPHA_TO[reg[j]];
+                }
+            if (q != 0)
+                continue;
+            /* store root and error location number indices */
+            root[count] = i;
+            loc[count] = k;
+            count++;
+        }
+        if (count != no_eras) {
+            printf("count = %d no_eras = %d\n lambda(x) is WRONG\n",count,no_eras);
+            count = -1;
+            goto finish;
+        }
 #if DEBUG >= 2
-    printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n");
-    for (i = 0; i < count; i++)
-      printf("%d ", loc[i]);
-    printf("\n");
+        printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n");
+        for (i = 0; i < count; i++)
+            printf("%d ", loc[i]);
+        printf("\n");
 #endif
 #endif
-  }
-  for(i=0;i<NROOTS+1;i++)
-    //    printf("%d  %d  %d\n",i,lambda[i],INDEX_OF[lambda[i]]);
-    b[i] = INDEX_OF[lambda[i]];
-  
-  /*
-   * Begin Berlekamp-Massey algorithm to determine error+erasure
-   * locator polynomial
-   */
-  r = no_eras;
-  el = no_eras;
-  while (++r <= NROOTS) {	/* r is the step number */
-    /* Compute discrepancy at the r-th step in poly-form */
-    discr_r = 0;
-    for (i = 0; i < r; i++){
-      if ((lambda[i] != 0) && (s[r-i-1] != A0)) {
-	discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r-i-1])];
-      }
     }
-    discr_r = INDEX_OF[discr_r];	/* Index form */
-    if (discr_r == A0) {
-      /* 2 lines below: B(x) <-- x*B(x) */
-      memmove(&b[1],b,NROOTS*sizeof(b[0]));
-      b[0] = A0;
-    } else {
-      /* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
-      t[0] = lambda[0];
-      for (i = 0 ; i < NROOTS; i++) {
-	if(b[i] != A0)
-	  t[i+1] = lambda[i+1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
-	else
-	  t[i+1] = lambda[i+1];
-      }
-      if (2 * el <= r + no_eras - 1) {
-	el = r + no_eras - el;
-	/*
-	 * 2 lines below: B(x) <-- inv(discr_r) *
-	 * lambda(x)
-	 */
-	for (i = 0; i <= NROOTS; i++)
-	  b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
-      } else {
-	/* 2 lines below: B(x) <-- x*B(x) */
-	memmove(&b[1],b,NROOTS*sizeof(b[0]));
-	b[0] = A0;
-      }
-      memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
-    }
-  }
-
-  /* Convert lambda to index form and compute deg(lambda(x)) */
-  deg_lambda = 0;
-  for(i=0;i<NROOTS+1;i++){
-    lambda[i] = INDEX_OF[lambda[i]];
-    if(lambda[i] != A0)
-      deg_lambda = i;
-  }
-  /* Find roots of the error+erasure locator polynomial by Chien search */
-  memcpy(&reg[1],&lambda[1],NROOTS*sizeof(reg[0]));
-  count = 0;		/* Number of roots of lambda(x) */
-  for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
-    q = 1; /* lambda[0] is always 0 */
-    for (j = deg_lambda; j > 0; j--){
-      if (reg[j] != A0) {
-	reg[j] = MODNN(reg[j] + j);
-	q ^= ALPHA_TO[reg[j]];
-      }
-    }
-    if (q != 0)
-      continue; /* Not a root */
-    /* store root (index-form) and error location number */
-#if DEBUG>=2
-    printf("count %d root %d loc %d\n",count,i,k);
-#endif
-    root[count] = i;
-    loc[count] = k;
-    /* If we've already found max possible roots,
-     * abort the search to save time
-     */
-    if(++count == deg_lambda)
-      break;
-  }
-  if (deg_lambda != count) {
-    /*
-     * deg(lambda) unequal to number of roots => uncorrectable
-     * error detected
-     */
-    count = -1;
-    goto finish;
-  }
-  /*
-   * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
-   * x**NROOTS). in index form. Also find deg(omega).
-   */
-  deg_omega = deg_lambda-1;
-  for (i = 0; i <= deg_omega;i++){
-    tmp = 0;
-    for(j=i;j >= 0; j--){
-      if ((s[i - j] != A0) && (lambda[j] != A0))
-	tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
-    }
-    omega[i] = INDEX_OF[tmp];
-  }
-
-  /*
-   * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
-   * inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
-   */
-  for (j = count-1; j >=0; j--) {
-    num1 = 0;
-    for (i = deg_omega; i >= 0; i--) {
-      if (omega[i] != A0)
-	num1  ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
-    }
-    num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
-    den = 0;
+    for(i=0;i<NROOTS+1;i++)
+        b[i] = INDEX_OF[lambda[i]];
     
-    /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
-    for (i = min(deg_lambda,NROOTS-1) & ~1; i >= 0; i -=2) {
-      if(lambda[i+1] != A0)
-	den ^= ALPHA_TO[MODNN(lambda[i+1] + i * root[j])];
+    /*
+     * Begin Berlekamp-Massey algorithm to determine error+erasure
+     * locator polynomial
+     */
+    r = no_eras;
+    el = no_eras;
+    while (++r <= NROOTS) {	/* r is the step number */
+        /* Compute discrepancy at the r-th step in poly-form */
+        discr_r = 0;
+        for (i = 0; i < r; i++){
+            if ((lambda[i] != 0) && (s[r-i-1] != A0)) {
+                discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r-i-1])];
+            }
+        }
+        discr_r = INDEX_OF[discr_r];	/* Index form */
+        if (discr_r == A0) {
+            /* 2 lines below: B(x) <-- x*B(x) */
+            memmove(&b[1],b,NROOTS*sizeof(b[0]));
+            b[0] = A0;
+        } else {
+            /* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
+            t[0] = lambda[0];
+            for (i = 0 ; i < NROOTS; i++) {
+                if(b[i] != A0)
+                    t[i+1] = lambda[i+1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
+                else
+                    t[i+1] = lambda[i+1];
+            }
+            if (2 * el <= r + no_eras - 1) {
+                el = r + no_eras - el;
+                /*
+                 * 2 lines below: B(x) <-- inv(discr_r) *
+                 * lambda(x)
+                 */
+                for (i = 0; i <= NROOTS; i++)
+                    b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
+            } else {
+                /* 2 lines below: B(x) <-- x*B(x) */
+                memmove(&b[1],b,NROOTS*sizeof(b[0]));
+                b[0] = A0;
+            }
+            memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
+        }
     }
-#if DEBUG >= 1
-    if (den == 0) {
-      printf("\n ERROR: denominator = 0\n");
-      count = -1;
-      goto finish;
+    
+    /* Convert lambda to index form and compute deg(lambda(x)) */
+    deg_lambda = 0;
+    for(i=0;i<NROOTS+1;i++){
+        lambda[i] = INDEX_OF[lambda[i]];
+        if(lambda[i] != A0)
+            deg_lambda = i;
     }
+    /* Find roots of the error+erasure locator polynomial by Chien search */
+    memcpy(&reg[1],&lambda[1],NROOTS*sizeof(reg[0]));
+    count = 0;		/* Number of roots of lambda(x) */
+    for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
+        q = 1; /* lambda[0] is always 0 */
+        for (j = deg_lambda; j > 0; j--){
+            if (reg[j] != A0) {
+                reg[j] = MODNN(reg[j] + j);
+                q ^= ALPHA_TO[reg[j]];
+            }
+        }
+        if (q != 0)
+            continue; /* Not a root */
+        /* store root (index-form) and error location number */
+#if DEBUG>=2
+        printf("count %d root %d loc %d\n",count,i,k);
 #endif
-    /* Apply error to data */
-    if (num1 != 0 && loc[j] >= PAD) {
-      data[loc[j]-PAD] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
+        root[count] = i;
+        loc[count] = k;
+        /* If we've already found max possible roots,
+         * abort the search to save time
+         */
+        if(++count == deg_lambda)
+            break;
     }
-  }
- finish:
-  if(eras_pos != NULL){
-    for(i=0;i<count;i++)
-      eras_pos[i] = loc[i];
-  }
-  return count;
+    if (deg_lambda != count) {
+        /*
+         * deg(lambda) unequal to number of roots => uncorrectable
+         * error detected
+         */
+        count = -1;
+        goto finish;
+    }
+    /*
+     * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
+     * x**NROOTS). in index form. Also find deg(omega).
+     */
+    deg_omega = 0;
+    for (i = 0; i < NROOTS;i++){
+        tmp = 0;
+        j = (deg_lambda < i) ? deg_lambda : i;
+        for(;j >= 0; j--){
+            if ((s[i - j] != A0) && (lambda[j] != A0))
+                tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
+        }
+        if(tmp != 0)
+            deg_omega = i;
+        omega[i] = INDEX_OF[tmp];
+    }
+    omega[NROOTS] = A0;
+    
+    /*
+     * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
+     * inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
+     */
+    for (j = count-1; j >=0; j--) {
+        num1 = 0;
+        for (i = deg_omega; i >= 0; i--) {
+            if (omega[i] != A0)
+                num1  ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
+        }
+        num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
+        den = 0;
+        
+        /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
+        for (i = min(deg_lambda,NROOTS-1) & ~1; i >= 0; i -=2) {
+            if(lambda[i+1] != A0)
+                den ^= ALPHA_TO[MODNN(lambda[i+1] + i * root[j])];
+        }
+        if (den == 0) {
+#if DEBUG >= 1
+            printf("\n ERROR: denominator = 0\n");
+#endif
+            count = -1;
+            goto finish;
+        }
+        /* Apply error to data */
+        if (num1 != 0) {
+            data[loc[j]] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
+        }
+    }
+finish:
+    if(eras_pos != NULL){
+        for(i=0;i<count;i++)
+            eras_pos[i] = loc[i];
+    }
+    return count;
 }
diff --git a/libm65/demod64a.f90 b/libm65/demod64a.f90
index 0f37095a5..3b8ee03ce 100644
--- a/libm65/demod64a.f90
+++ b/libm65/demod64a.f90
@@ -33,10 +33,12 @@ subroutine demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
   do j=1,63
      s1=-1.e30
      fsum=0.
+     psum=0.
      do i=1,64
         x=min(afac*s3(i,j)/ave,50.d0)
         fs(i)=exp(x)
         fsum=fsum+fs(i)
+        psum=psum + s3(i,j)
         if(s3(i,j).gt.s1) then
            s1=s3(i,j)
            i1=i                              !Most reliable
@@ -50,8 +52,10 @@ subroutine demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
            i2=i                              !Second most reliable
         endif
      enddo
-     p1=fs(i1)/fsum                          !Normalized probabilities
-     p2=fs(i2)/fsum
+!     p1=fs(i1)/fsum                          !Normalized probabilities
+!     p2=fs(i2)/fsum
+     p1=s1/psum
+     p2=s2/psum
      mrsym(j)=i1-1
      mr2sym(j)=i2-1
      mrprob(j)=scale*p1
@@ -66,7 +70,8 @@ subroutine demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
      sum=sum+mrprob(j)
      if(mrprob(j).le.5) nlow=nlow+1
   enddo
-  ntest=sum/63
+!  ntest=sum/63
+  ntest=sum
 
   return
 end subroutine demod64a
diff --git a/libm65/encode_rs.c b/libm65/encode_rs.c
index 911bc36dd..9d56d0bf1 100644
--- a/libm65/encode_rs.c
+++ b/libm65/encode_rs.c
@@ -13,24 +13,19 @@
 #endif
 
 void ENCODE_RS(
-#ifdef FIXED
-DTYPE *data, DTYPE *bb,int pad){
-#else
-void *p,DTYPE *data, DTYPE *bb){
+#ifndef FIXED
+void *p,
+#endif
+DTYPE *data, DTYPE *bb){
+#ifndef FIXED
   struct rs *rs = (struct rs *)p;
 #endif
   int i, j;
   DTYPE feedback;
 
-#ifdef FIXED
-  /* Check pad parameter for validity */
-  if(pad < 0 || pad >= NN)
-    return;
-#endif
-
   memset(bb,0,NROOTS*sizeof(DTYPE));
 
-  for(i=0;i<NN-NROOTS-PAD;i++){
+  for(i=0;i<NN-NROOTS;i++){
     feedback = INDEX_OF[data[i] ^ bb[0]];
     if(feedback != A0){      /* feedback term is non-zero */
 #ifdef UNNORMALIZED
diff --git a/libm65/extract.f90 b/libm65/extract.f90
index 6e3457ee0..f59357ec1 100644
--- a/libm65/extract.f90
+++ b/libm65/extract.f90
@@ -6,26 +6,36 @@ subroutine extract(s3,nadd,ncount,nhist,decoded,ltext)
   integer era(51),dat4(12),indx(64)
   integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
   logical first,ltext
+  integer correct(63),itmp(63)
+  integer param(0:8)
+  integer h0(0:11),d0(0:11),ne(0:11)
+  real r0(0:11)
+  common/test001/s3a(64,63),mrs(63),mrs2(63)        !### TEST ONLY ###
+
+!          0  1  2  3  4  5  6  7  8  9 10 11
+  data h0/41,42,43,43,44,45,46,47,48,48,49,49/
+  data d0/71,72,73,74,76,77,78,80,81,82,83,83/
+!            0    1    2    3    4    5    6    7    8    9   10   11
+  data r0/0.70,0.72,0.74,0.76,0.78,0.80,0.82,0.84,0.86,0.88,0.90,0.90/
+
   data first/.true./,nsec1/0/
   save
 
   nfail=0
-1 continue
-!      call timer('demod64a',0)
-  call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
-!      call timer('demod64a',1)
-  if(ntest.lt.50 .or. nlow.gt.20) then
-     ncount=-999                         !Flag bad data
-     go to 900
-  endif
+  call pctile(s3,tmp,4032,50,base)     ! ### or, use ave from demod64a
+  s3=s3/base
+  s3a=s3
+1 call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
+!  if(ntest.lt.50 .or. nlow.gt.20) then
+!     ncount=-999                         !Flag bad data
+!     go to 900
+!  endif
   call chkhist(mrsym,nhist,ipk)
 
   if(nhist.ge.20) then
      nfail=nfail+1
      call pctile(s3,tmp,4032,50,base)     ! ### or, use ave from demod64a
-     do j=1,63
-        s3(ipk,j)=base
-     enddo
+     s3(ipk,1:63)=base
      if(nfail.gt.30) then
         decoded='                      '
         ncount=-1
@@ -34,74 +44,92 @@ subroutine extract(s3,nadd,ncount,nhist,decoded,ltext)
      go to 1
   endif
 
+  mrs=mrsym
+  mrs2=mr2sym
+
   call graycode(mrsym,63,-1)
   call interleave63(mrsym,-1)
   call interleave63(mrprob,-1)
 
-  ndec=1
-  nemax=30                              !Was 200 (30)
-  maxe=8
-  xlambda=13.0                          !Was 12
+  call graycode(mr2sym,63,-1)
+  call interleave63(mr2sym,-1)
+  call interleave63(mr2prob,-1)
 
-  if(ndec.eq.1) then
-     call graycode(mr2sym,63,-1)
-     call interleave63(mr2sym,-1)
-     call interleave63(mr2prob,-1)
+  ntrials=10000
+  naggressive=10
 
-     nsec1=nsec1+1
-     write(22,rec=1) nsec1,xlambda,maxe,200,mrsym,mrprob,mr2sym,mr2prob
-     call flush(22)
-!         call timer('kvasd   ',0)
-!#ifdef UNIX
-!         iret=system('./kvasd -q > dev_null')
-!#else
-     iret=system('kvasd -q > dev_null')
-!#endif
-!         call timer('kvasd   ',1)
-     if(iret.ne.0) then
-        if(first) write(*,1000) iret
-1000    format('Error in KV decoder, or no KV decoder present.'/     &
-               'Return code:',i8,'.  Will use BM algorithm.')
-        ndec=0
-        first=.false.
-        go to 20
-     endif
+  ntry=0
+  param=0
 
-     read(22,rec=2) nsec2,ncount,dat4
-     j=nsec2                !Silence compiler warning
-     decoded='                      '
-     ltext=.false.
-     if(ncount.ge.0) then
-        call unpackmsg(dat4,decoded)     !Unpack the user message
-        if(iand(dat4(10),8).ne.0) ltext=.true.
-        do i=2,12
-           if(dat4(i).ne.dat4(1)) go to 20
-        enddo
-        write(13,*) 'Bad decode?',nhist,nfail,ipk,' ',dat4,decoded
-        ncount=-1           !Suppress supposedly bogus decodes
-        decoded='                      '
-     endif
-  endif
-20 if(ndec.eq.0) then
-     call indexx(63,mrprob,indx)
-     do i=1,nemax
-        j=indx(i)
-        if(mrprob(j).gt.120) then
-           ne2=i-1
-           go to 2
-        endif
-        era(i)=j-1
-     enddo
-     ne2=nemax
-2    decoded='                      '
-     do nerase=0,ne2,2
-        call rs_decode(mrsym,era,nerase,dat4,ncount)
-        if(ncount.ge.0) then
-           call unpackmsg(dat4,decoded)
-           go to 900
-        endif
-     enddo
+  call timer('ftrsd   ',0)
+  call ftrsd2(mrsym,mrprob,mr2sym,mr2prob,ntrials,correct,param,ntry)
+  call timer('ftrsd   ',1)
+  ncandidates=param(0)
+  nhard=param(1)
+  nsoft=param(2)
+  nerased=param(3)
+  rtt=0.001*param(4)
+  ntotal=param(5)
+  qual=0.001*param(7)
+  nd0=81
+  r00=0.87
+  if(naggressive.eq.10) then
+     nd0=83
+     r00=0.90
   endif
+  if(ntotal.le.nd0 .and. rtt.le.r00) nft=1
+  n=naggressive
+  if(nhard.gt.50) nft=0
+  if(nhard.gt.h0(n)) nft=0
+  if(ntotal.gt.d0(n)) nft=0
+  if(rtt.gt.r0(n)) nft=0
 
-900 return
+  ncount=-1
+  decoded='                      '
+  ltext=.false.
+  if(nft.gt.0) then
+! Turn the corrected symbol array into channel symbols for subtraction;
+! pass it back to jt65a via common block "chansyms65".
+     do i=1,12
+        dat4(i)=correct(13-i)
+     enddo
+     do i=1,63
+       itmp(i)=correct(64-i)
+     enddo
+     correct(1:63)=itmp(1:63)
+     call interleave63(correct,63,1)
+     call graycode65(correct,63,1)
+     call unpackmsg(dat4,decoded)     !Unpack the user message
+     ncount=0
+     if(iand(dat4(10),8).ne.0) ltext=.true.
+  endif
+900 continue
+  if(nft.eq.1 .and. nhard.lt.0) decoded='                      '
+!  write(81,3001) naggressive,ncandidates,nhard,ntotal,rtt,qual,decoded
+!3001 format(i2,i6,i3,i4,2f8.2,2x,a22)
+
+  return
 end subroutine extract
+
+subroutine getpp(workdat,p)
+
+  integer workdat(63)
+  integer a(63)
+  common/test001/s3a(64,63),mrs(63),mrs2(63)
+
+  a(1:63)=workdat(63:1:-1)
+  call interleave63(a,1)
+  call graycode(a,63,1,a)
+
+  psum=0.
+  do j=1,63
+     i=a(j)+1
+     x=s3a(i,j)
+     s3a(i,j)=0.
+     psum=psum + x
+     s3a(i,j)=x
+  enddo
+  p=psum/63.0
+
+  return
+end subroutine getpp
diff --git a/libm65/ftrsd2.c b/libm65/ftrsd2.c
new file mode 100644
index 000000000..241312483
--- /dev/null
+++ b/libm65/ftrsd2.c
@@ -0,0 +1,213 @@
+/*
+ ftrsd2.c
+ 
+ A soft-decision decoder for the JT65 (63,12) Reed-Solomon code.
+ 
+ This decoding scheme is built around Phil Karn's Berlekamp-Massey
+ errors and erasures decoder. The approach is inspired by a number of
+ publications, including the stochastic Chase decoder described
+ in "Stochastic Chase Decoding of Reed-Solomon Codes", by Leroux et al.,
+ IEEE Communications Letters, Vol. 14, No. 9, September 2010 and
+ "Soft-Decision Decoding of Reed-Solomon Codes Using Successive Error-
+ and-Erasure Decoding," by Soo-Woong Lee and B. V. K. Vijaya Kumar.
+ 
+ Steve Franke K9AN and Joe Taylor K1JT
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <time.h>
+#include <string.h>
+#include "rs2.h"
+
+static void *rs;
+void getpp_(int workdat[], float *pp);
+
+void ftrsd2_(int mrsym[], int mrprob[], int mr2sym[], int mr2prob[], 
+	     int* ntrials0, int correct[], int param[], int ntry[])
+{
+  int rxdat[63], rxprob[63], rxdat2[63], rxprob2[63];
+  int workdat[63];
+  int indexes[63];
+  int era_pos[51];
+  int i, j, numera, nerr, nn=63;
+  int ntrials = *ntrials0;
+  int nhard=0,nhard_min=32768,nsoft=0,nsoft_min=32768;
+  int ntotal=0,ntotal_min=32768,ncandidates;
+  int nera_best=0;
+  float pp,pp1,pp2;
+  static unsigned int nseed;
+
+// Power-percentage symbol metrics - composite gnnf/hf 
+  int perr[8][8] = {
+    { 4,      9,     11,     13,     14,     14,     15,     15},
+    { 2,     20,     20,     30,     40,     50,     50,     50},
+    { 7,     24,     27,     40,     50,     50,     50,     50},
+    {13,     25,     35,     46,     52,     70,     50,     50},
+    {17,     30,     42,     54,     55,     64,     71,     70},
+    {25,     39,     48,     57,     64,     66,     77,     77},
+    {32,     45,     54,     63,     66,     75,     78,     83},
+    {51,     58,     57,     66,     72,     77,     82,     86}};
+
+    
+// Initialize the KA9Q Reed-Solomon encoder/decoder
+  unsigned int symsize=6, gfpoly=0x43, fcr=3, prim=1, nroots=51;
+  rs=init_rs_int(symsize, gfpoly, fcr, prim, nroots, 0);
+
+// Reverse the received symbol vectors for BM decoder
+  for (i=0; i<63; i++) {
+    rxdat[i]=mrsym[62-i];
+    rxprob[i]=mrprob[62-i];
+    rxdat2[i]=mr2sym[62-i];
+    rxprob2[i]=mr2prob[62-i];
+  }
+    
+// Sort rxprob to find indexes of the least reliable symbols
+  int k, pass, tmp, nsym=63;
+  int probs[63];
+  for (i=0; i<63; i++) {
+    indexes[i]=i;
+    probs[i]=rxprob[i];
+  }
+  for (pass = 1; pass <= nsym-1; pass++) {
+    for (k = 0; k < nsym - pass; k++) {
+      if( probs[k] < probs[k+1] ) {
+        tmp = probs[k];
+        probs[k] = probs[k+1];
+        probs[k+1] = tmp;
+        tmp = indexes[k];
+        indexes[k] = indexes[k+1];
+        indexes[k+1] = tmp;
+      }
+    }
+  }
+  
+// See if we can decode using BM HDD, and calculate the syndrome vector.
+  memset(era_pos,0,51*sizeof(int));
+  numera=0;
+  memcpy(workdat,rxdat,sizeof(rxdat));
+  nerr=decode_rs_int(rs,workdat,era_pos,numera,1);
+  if( nerr >= 0 ) {
+    // Hard-decision decoding succeeded.  Save codeword and some parameters.
+    nhard=0;
+    for (i=0; i<63; i++) {
+      if( workdat[i] != rxdat[i] ) nhard=nhard+1;
+    }
+    memcpy(correct,workdat,63*sizeof(int));
+    param[0]=0;
+    param[1]=nhard;
+    param[2]=0;
+    param[3]=0;
+    param[4]=0;
+    param[5]=0;
+    param[7]=1000*1000;
+    ntry[0]=0;
+    return;
+  }
+
+/*
+Hard-decision decoding failed.  Try the FT soft-decision method.
+Generate random erasure-locator vectors and see if any of them
+decode. This will generate a list of "candidate" codewords.  The
+soft distance between each candidate codeword and the received 
+word is estimated by finding the largest (pp1) and second-largest 
+(pp2) outputs from a synchronized filter-bank operating on the 
+symbol spectra, and using these to decide which candidate 
+codeword is "best".
+*/
+
+  nseed=1;                                 //Seed for random numbers
+  float ratio;
+  int thresh, nsum;
+  int thresh0[63];
+  ncandidates=0;
+  nsum=0;
+  int ii,jj;
+  for (i=0; i<nn; i++) {
+    nsum=nsum+rxprob[i];
+    j = indexes[62-i];
+    ratio = (float)rxprob2[j]/((float)rxprob[j]+0.01);
+    ii = 7.999*ratio;
+    jj = (62-i)/8;
+    thresh0[i] = 1.3*perr[ii][jj];
+  }
+  if(nsum<=0) return;
+
+  pp1=0.0;
+  pp2=0.0;
+  for (k=1; k<=ntrials; k++) {
+    memset(era_pos,0,51*sizeof(int));
+    memcpy(workdat,rxdat,sizeof(rxdat));
+
+/* 
+Mark a subset of the symbols as erasures.
+Run through the ranked symbols, starting with the worst, i=0.
+NB: j is the symbol-vector index of the symbol with rank i.
+*/
+    numera=0;
+    for (i=0; i<nn; i++) {
+      j = indexes[62-i];
+      thresh=thresh0[i];
+      long int ir;
+
+// Generate a random number ir, 0 <= ir < 100 (see POSIX.1-2001 example).
+      nseed = nseed * 1103515245 + 12345;
+      ir = (unsigned)(nseed/65536) % 32768;
+      ir = (100*ir)/32768;
+
+      if((ir < thresh ) && numera < 51) {
+        era_pos[numera]=j;
+        numera=numera+1;
+      }
+    }
+
+    nerr=decode_rs_int(rs,workdat,era_pos,numera,0);        
+    if( nerr >= 0 ) {
+      // We have a candidate codeword.  Find its hard and soft distance from
+      // the received word.  Also find pp1 and pp2 from the full array 
+      // s3(64,63) of synchronized symbol spectra.
+      ncandidates=ncandidates+1;
+      nhard=0;
+      nsoft=0;
+      for (i=0; i<63; i++) {
+        if(workdat[i] != rxdat[i]) {
+          nhard=nhard+1;
+          if(workdat[i] != rxdat2[i]) {
+            nsoft=nsoft+rxprob[i];
+          }
+        }
+      }
+      nsoft=63*nsoft/nsum;
+      ntotal=nsoft+nhard;
+
+      getpp_(workdat,&pp);
+      if(pp>pp1) {
+        pp2=pp1;
+        pp1=pp;
+        nsoft_min=nsoft;
+        nhard_min=nhard;
+        ntotal_min=ntotal;
+        memcpy(correct,workdat,63*sizeof(int));
+        nera_best=numera;
+        ntry[0]=k;
+      } else {
+        if(pp>pp2 && pp!=pp1) pp2=pp;
+      }
+      if(nhard_min <= 41 && ntotal_min <= 71) break;
+    }
+    if(k == ntrials) ntry[0]=k;
+  }
+  
+  param[0]=ncandidates;
+  param[1]=nhard_min;
+  param[2]=nsoft_min;
+  param[3]=nera_best;
+  param[4]=1000.0*pp2/pp1;
+  param[5]=ntotal_min;
+  param[6]=ntry[0];
+  param[7]=1000.0*pp2;
+  param[8]=1000.0*pp1;
+  if(param[0]==0) param[2]=-1;
+  return;
+}
diff --git a/libm65/graycode65.f90 b/libm65/graycode65.f90
new file mode 100644
index 000000000..bb2c669ce
--- /dev/null
+++ b/libm65/graycode65.f90
@@ -0,0 +1,9 @@
+subroutine graycode65(dat,n,idir)
+
+  integer dat(n)
+  do i=1,n
+     dat(i)=igray(dat(i),idir)
+  enddo
+
+  return
+end subroutine graycode65
diff --git a/libm65/init_rs.c b/libm65/init_rs.c
index 7e1e09e5b..ce9e8cb70 100644
--- a/libm65/init_rs.c
+++ b/libm65/init_rs.c
@@ -30,30 +30,25 @@ void FREE_RS(void *p){
  * fcr = first root of RS code generator polynomial, index form
  * prim = primitive element to generate polynomial roots
  * nroots = RS code generator polynomial degree (number of roots)
- * pad = padding bytes at front of shortened block
  */
-void *INIT_RS(int symsize,int gfpoly,int fcr,int prim,
-	int nroots,int pad){
+void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned fcr,unsigned prim,
+		unsigned int nroots){
   struct rs *rs;
   int i, j, sr,root,iprim;
 
-  /* Check parameter ranges */
-  if(symsize < 0 || symsize > 8*sizeof(DTYPE))
+  if(symsize > 8*sizeof(DTYPE))
     return NULL; /* Need version with ints rather than chars */
 
-  if(fcr < 0 || fcr >= (1<<symsize))
+  if(fcr >= (1<<symsize))
     return NULL;
-  if(prim <= 0 || prim >= (1<<symsize))
+  if(prim == 0 || prim >= (1<<symsize))
     return NULL;
-  if(nroots < 0 || nroots >= (1<<symsize))
+  if(nroots >= (1<<symsize))
     return NULL; /* Can't have more roots than symbol values! */
-  if(pad < 0 || pad >= ((1<<symsize) -1 - nroots))
-    return NULL; /* Too much padding */
 
   rs = (struct rs *)calloc(1,sizeof(struct rs));
   rs->mm = symsize;
   rs->nn = (1<<symsize)-1;
-  rs->pad = pad;
 
   rs->alpha_to = (DTYPE *)malloc(sizeof(DTYPE)*(rs->nn+1));
   if(rs->alpha_to == NULL){
diff --git a/libm65/int.h b/libm65/int.h
index 056241e12..ada5bfd4c 100644
--- a/libm65/int.h
+++ b/libm65/int.h
@@ -7,19 +7,18 @@
 
 /* Reed-Solomon codec control block */
 struct rs {
-  int mm;              /* Bits per symbol */
-  int nn;              /* Symbols per block (= (1<<mm)-1) */
-  DTYPE *alpha_to;     /* log lookup table */
-  DTYPE *index_of;     /* Antilog lookup table */
-  DTYPE *genpoly;      /* Generator polynomial */
-  int nroots;     /* Number of generator roots = number of parity symbols */
-  int fcr;        /* First consecutive root, index form */
-  int prim;       /* Primitive element, index form */
-  int iprim;      /* prim-th root of 1, index form */
-  int pad;        /* Padding bytes in shortened block */
+  unsigned int mm;   /* Bits per symbol */
+  unsigned int nn;   /* Symbols per block (= (1<<mm)-1) */
+  int *alpha_to;      /* log lookup table */
+  int *index_of;      /* Antilog lookup table */
+  int *genpoly;       /* Generator polynomial */
+  unsigned int nroots;     /* Number of generator roots = number of parity symbols */
+  unsigned int fcr;        /* First consecutive root, index form */
+  unsigned int prim;       /* Primitive element, index form */
+  unsigned int iprim;      /* prim-th root of 1, index form */
 };
 
-static int modnn(struct rs *rs,int x){
+static inline int modnn(struct rs *rs,int x){
   while (x >= rs->nn) {
     x -= rs->nn;
     x = (x >> rs->mm) + (x & rs->nn);
@@ -33,12 +32,10 @@ static int modnn(struct rs *rs,int x){
 #define ALPHA_TO (rs->alpha_to) 
 #define INDEX_OF (rs->index_of)
 #define GENPOLY (rs->genpoly)
-//#define NROOTS (rs->nroots)
-#define NROOTS (51)
+#define NROOTS (rs->nroots)
 #define FCR (rs->fcr)
 #define PRIM (rs->prim)
 #define IPRIM (rs->iprim)
-#define PAD (rs->pad)
 #define A0 (NN)
 
 #define ENCODE_RS encode_rs_int
@@ -47,9 +44,9 @@ static int modnn(struct rs *rs,int x){
 #define FREE_RS free_rs_int
 
 void ENCODE_RS(void *p,DTYPE *data,DTYPE *parity);
-int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras);
-void *INIT_RS(int symsize,int gfpoly,int fcr,
-		   int prim,int nroots,int pad);
+int DECODE_RS(void *p,DTYPE *data,int *eras_pos,int no_eras, int calc_syn);
+void *INIT_RS(unsigned int symsize,unsigned int gfpoly,unsigned int fcr,
+		   unsigned int prim,unsigned int nroots);
 void FREE_RS(void *p);
 
 
diff --git a/libm65/rs2.h b/libm65/rs2.h
new file mode 100644
index 000000000..c2b807d15
--- /dev/null
+++ b/libm65/rs2.h
@@ -0,0 +1,16 @@
+/* User include file for the Reed-Solomon codec
+ * Copyright 2002, Phil Karn KA9Q
+ * May be used under the terms of the GNU General Public License (GPL)
+ */
+
+/* General purpose RS codec, integer symbols */
+void encode_rs_int(void *rs,int *data,int *parity);
+int decode_rs_int(void *rs,int *data,int *eras_pos,int no_eras, int calc_syn);
+void *init_rs_int(int symsize,int gfpoly,int fcr,
+		  int prim,int nroots,int pad);
+void free_rs_int(void *rs);
+
+/* Tables to map from conventional->dual (Taltab) and
+ * dual->conventional (Tal1tab) bases
+ */
+extern unsigned char Taltab[],Tal1tab[];
diff --git a/mainwindow.cpp b/mainwindow.cpp
index 19a5a3707..c11b80783 100644
--- a/mainwindow.cpp
+++ b/mainwindow.cpp
@@ -33,7 +33,7 @@ TxTune*    g_pTxTune = NULL;
 QSharedMemory mem_m65("mem_m65");
 
 QString rev="$Rev$";                            //Must update by hand ????
-QString Program_Title_Version="  MAP65   v2.5, r" + rev.mid(6,4) +
+QString Program_Title_Version="  MAP65   v2.6, r" + rev.mid(6,4) +
                               "    by K1JT";
 
 extern const int RxDataFrequency = 96000;