From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Fri, 29 May 2015 06:33:00 +0000 (+0000)
Subject: prog cmd line testing, can now compare between SSB and cohpsk
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=75e27716165f1054791b16b96fa78d5ccbd0e0cf;p=freetel-svn-tracking.git

prog cmd line testing, can now compare between SSB and cohpsk

git-svn-id: https://svn.code.sf.net/p/freetel/code@2163 01035d8c-6547-0410-b346-abe4f91aad63
---

diff --git a/codec2-dev/octave/cellmodem.m b/codec2-dev/octave/cellmodem.m
index 5bf7626f..7477d06f 100644
--- a/codec2-dev/octave/cellmodem.m
+++ b/codec2-dev/octave/cellmodem.m
@@ -6,12 +6,18 @@
 % Ideas:
 %   + insert low rate codec
 %   + generate bunch of symbols, run through codec, measure MSE, choose best set
+%   + measure probablility of error, "distance" from other symbols
+%   + can we use VQ training algorithm for this?  Start with 2 symbols,
+%     pass through channel, measure MSE, split again?
 %   + start with cmd line version of codec, frame synchronous
 %   + add symbol timing estimator later
 %   + try different frame rates
 %   + simulate impairments like HP/LP filtering.  Can we correct for this?
 %   + pilots symbols so we can use energy as well?
 %   + set of F0 as well
+%   + LSP quantisers preferrentially preserve peaks, so three peaks is a
+%     reasonable signal set.  Modulate position and bandwidth to creat 
+%     symbol set.  Maybe 50 or 100 Hz grid for LSPs, evaluate that some how.
 
 graphics_toolkit ("gnuplot");
 lsp;
@@ -36,6 +42,8 @@ for m=1:l
 end
 s = filter(1, a, ex);
 
+% insert codec here (e.g. write to file, run codec, read file)
+
 % extract received symbol from channel
 
 a_ = lpcauto(s, lpc_order)
diff --git a/codec2-dev/octave/make_ssbfilt.m b/codec2-dev/octave/make_ssbfilt.m
new file mode 100644
index 00000000..b1bb3474
--- /dev/null
+++ b/codec2-dev/octave/make_ssbfilt.m
@@ -0,0 +1,36 @@
+% make_ssbfilt.m
+% David Rowe May 2015
+%
+% creates SSB filter FIR coeffs
+
+graphics_toolkit ("gnuplot");
+
+ssbfilt_n = 100;
+Fs = 8000;
+
+ssbfilt_coeff  = fir2(ssbfilt_n,[0 400 600 2200 2600 4000]/(Fs/2),[0.001 0.001 1 1 0.001 0.001]);
+
+figure(1)
+clf;
+h = freqz(ssbfilt_coeff,1,Fs/2);
+plot(20*log10(abs(h)))
+grid minor
+
+% save coeffs to a C header file
+
+f=fopen("../src/ssbfilt_coeff.h","wt");
+fprintf(f,"/* 600 - 2600 Hz FIR filter coeffs */\n");
+fprintf(f,"/* Generated by make_ssbfilt Octave script */\n");
+
+fprintf(f,"\n#define SSBFILT_N %d\n\n", ssbfilt_n);
+
+fprintf(f,"float ssbfilt_coeff[]={\n");
+for r=1:ssbfilt_n
+  if r < ssbfilt_n
+    fprintf(f, "  %f,\n",  ssbfilt_coeff(r));
+  else
+    fprintf(f, "  %f\n};", ssbfilt_coeff(r));
+  end
+end
+
+fclose(f);
diff --git a/codec2-dev/octave/test_cohpsk_ch.m b/codec2-dev/octave/test_cohpsk_ch.m
index 73446641..19634ddf 100644
--- a/codec2-dev/octave/test_cohpsk_ch.m
+++ b/codec2-dev/octave/test_cohpsk_ch.m
@@ -5,7 +5,7 @@
 
 Nc=7; Nd=2;
 
-load ../build_linux/unittest/test_cohpsk_ch_out.txt
+load ../build_linux/src/test_cohpsk_ch_out.txt
   
 figure(3)
 clf;
diff --git a/codec2-dev/raw/ve9qrp_10s.raw b/codec2-dev/raw/ve9qrp_10s.raw
new file mode 100644
index 00000000..0d3bd3b9
Binary files /dev/null and b/codec2-dev/raw/ve9qrp_10s.raw differ
diff --git a/codec2-dev/src/CMakeLists.txt b/codec2-dev/src/CMakeLists.txt
index fbf7e3e6..05ce4675 100644
--- a/codec2-dev/src/CMakeLists.txt
+++ b/codec2-dev/src/CMakeLists.txt
@@ -278,7 +278,7 @@ target_link_libraries(fm_demod ${CMAKE_REQUIRED_LIBRARIES})
 add_executable(cohpsk_mod cohpsk_mod.c)
 target_link_libraries(cohpsk_mod ${CMAKE_REQUIRED_LIBRARIES} codec2)
 
-add_executable(cohpsk_demod cohpsk_demod.c)
+add_executable(cohpsk_demod cohpsk_demod.c octave.c)
 target_link_libraries(cohpsk_demod ${CMAKE_REQUIRED_LIBRARIES} codec2)
 
 add_executable(cohpsk_get_test_bits cohpsk_get_test_bits.c)
diff --git a/codec2-dev/src/c2dec.c b/codec2-dev/src/c2dec.c
index 880db0b3..2330c910 100644
--- a/codec2-dev/src/c2dec.c
+++ b/codec2-dev/src/c2dec.c
@@ -51,12 +51,13 @@ int main(int argc, char *argv[])
     FILE          *fber = NULL;
     short         *buf;
     unsigned char *bits, *prev_bits;
+    int           *unpacked_bits;
     int            nsam, nbit, nbyte, i, byte, frames, bits_proc, bit_errors, error_mode;
     int            nstart_bit, nend_bit, bit_rate;
     int            state, next_state;
     float          ber, r, burst_length, burst_period, burst_timer, ber_est;
     unsigned char  mask;
-    int            natural, dump;
+    int            natural, dump, unpacked, bit, ret;
 
     char* opt_string = "h:";
     struct option long_options[] = {
@@ -65,6 +66,7 @@ int main(int argc, char *argv[])
         { "endbit", required_argument, NULL, 0 },
         { "berfile", required_argument, NULL, 0 },
         { "natural", no_argument, &natural, 1 },
+        { "unpacked", no_argument, &unpacked, 1 },
         #ifdef DUMP
         { "dump", required_argument, &dump, 1 },
         #endif
@@ -88,10 +90,10 @@ int main(int argc, char *argv[])
 	mode = CODEC2_MODE_1300;
     else if (strcmp(argv[1],"1200") == 0)
 	mode = CODEC2_MODE_1200;
-    else if (strcmp(argv[1],"650") == 0)
-	mode = CODEC2_MODE_650;
+    else if (strcmp(argv[1],"700") == 0)
+	mode = CODEC2_MODE_700;
     else {
-	fprintf(stderr, "Error in mode: %s.  Must be 3200, 2400, 1600, 1400, 1300, 1200, or 650\n", argv[1]);
+	fprintf(stderr, "Error in mode: %s.  Must be 3200, 2400, 1600, 1400, 1300, 1200, or 700\n", argv[1]);
 	exit(1);
     }
     bit_rate = atoi(argv[1]);
@@ -122,6 +124,7 @@ int main(int argc, char *argv[])
     buf = (short*)malloc(nsam*sizeof(short));
     nbyte = (nbit + 7) / 8;
     bits = (unsigned char*)malloc(nbyte*sizeof(char));
+    unpacked_bits = (int*)malloc(nbit*sizeof(int));
     prev_bits = (unsigned char*)malloc(nbyte*sizeof(char));
     frames = bit_errors = bits_proc = 0;
     nstart_bit = 0;
@@ -172,12 +175,18 @@ int main(int argc, char *argv[])
     codec2_set_natural_or_gray(codec2, !natural);
     //printf("%d %d\n", nstart_bit, nend_bit);
  
-    while(fread(bits, sizeof(char), nbyte, fin) == (size_t)nbyte) {
+    if (unpacked)
+        ret = (fread(unpacked_bits, sizeof(int), nbit, fin) == (size_t)nbit);
+    else
+        ret = (fread(bits, sizeof(char), nbyte, fin) == (size_t)nbyte);
+      
+    while(ret) {
 	frames++;
 
-        // apply bit errors, MSB of byte 0 is bit 0 in frame */
+        // apply bit errors, MSB of byte 0 is bit 0 in frame, only works in packed mode
         
 	if ((error_mode == UNIFORM) || (error_mode == UNIFORM_RANGE)) {
+            assert(unpacked == 0);
 	    for(i=nstart_bit; i<nend_bit+1; i++) {
 		r = (float)rand()/RAND_MAX;
 		if (r < ber) {
@@ -193,6 +202,7 @@ int main(int argc, char *argv[])
 	}
 
 	if (error_mode == TWO_STATE) {
+            assert(unpacked == 0);
             burst_timer += (float)nbit/bit_rate;
             fprintf(stderr, "burst_timer: %f  state: %d\n", burst_timer, state);
 
@@ -241,14 +251,36 @@ int main(int argc, char *argv[])
         else
             ber_est = 0.0;
             
+        if (unpacked) {
+            /* pack bits, MSB received first  */
+
+            bit = 7; byte = 0;
+            memset(bits, 0, nbyte);
+            for(i=0; i<nbit; i++) {
+                bits[byte] |= (unpacked_bits[i] << bit);
+                bit--;
+                if (bit < 0) {
+                    bit = 7;
+                    byte++;
+                }
+            }
+       }
+
 	codec2_decode_ber(codec2, buf, bits, ber_est);
  	fwrite(buf, sizeof(short), nsam, fout);
+
 	//if this is in a pipeline, we probably don't want the usual
         //buffering to occur
+
         if (fout == stdout) fflush(stdout);
         if (fin == stdin) fflush(stdin);         
 
         memcpy(prev_bits, bits, nbyte);
+
+        if (unpacked)
+            ret = (fread(unpacked_bits, sizeof(int), nbit, fin) == (size_t)nbit);
+        else
+            ret = (fread(bits, sizeof(char), nbyte, fin) == (size_t)nbyte);
     }
 
     if (error_mode)
@@ -258,6 +290,7 @@ int main(int argc, char *argv[])
 
     free(buf);
     free(bits);
+    free(unpacked_bits);
     fclose(fin);
     fclose(fout);
 
diff --git a/codec2-dev/src/c2enc.c b/codec2-dev/src/c2enc.c
index 1aa57563..5f5fddab 100644
--- a/codec2-dev/src/c2enc.c
+++ b/codec2-dev/src/c2enc.c
@@ -41,10 +41,12 @@ int main(int argc, char *argv[])
     FILE          *fout;
     short         *buf;
     unsigned char *bits;
-    int            nsam, nbit, nbyte, gray;
- 
+    int            nsam, nbit, nbyte, gray, unpacked;
+    int           *unpacked_bits;
+    int            bit, byte,i;
+
     if (argc < 4) {
-	printf("usage: c2enc 3200|2400|1600|1400|1300|1200|650 InputRawspeechFile OutputBitFile [--natural]\n");
+	printf("usage: c2enc 3200|2400|1600|1400|1300|1200|700 InputRawspeechFile OutputBitFile [--natural] [--unpacked]\n");
 	printf("e.g    c2enc 1400 ../raw/hts1a.raw hts1a.c2\n");
 	printf("e.g    c2enc 1300 ../raw/hts1a.raw hts1a.c2 --natural\n");
 	exit(1);
@@ -62,10 +64,10 @@ int main(int argc, char *argv[])
 	mode = CODEC2_MODE_1300;
     else if (strcmp(argv[1],"1200") == 0)
 	mode = CODEC2_MODE_1200;
-    else if (strcmp(argv[1],"650") == 0)
-	mode = CODEC2_MODE_650;
+    else if (strcmp(argv[1],"700") == 0)
+	mode = CODEC2_MODE_700;
     else {
-	fprintf(stderr, "Error in mode: %s.  Must be 3200, 2400, 1600, 1400, 1300, 1200 or 650\n", argv[1]);
+	fprintf(stderr, "Error in mode: %s.  Must be 3200, 2400, 1600, 1400, 1300, 1200 or 700\n", argv[1]);
 	exit(1);
     }
 
@@ -90,6 +92,7 @@ int main(int argc, char *argv[])
     nbyte = (nbit + 7) / 8;
 
     bits = (unsigned char*)malloc(nbyte*sizeof(char));
+    unpacked_bits = (int*)malloc(nbit*sizeof(int));
     
     if (argc == 5) {
         if (strcmp(argv[4], "--natural") == 0)
@@ -97,13 +100,37 @@ int main(int argc, char *argv[])
         else
             gray = 1;
         codec2_set_natural_or_gray(codec2, gray);
+        if (strcmp(argv[4], "--unpacked") == 0) {
+            unpacked = 1;
+        }
+        else
+            unpacked = 0;
     }
 
     while(fread(buf, sizeof(short), nsam, fin) == (size_t)nsam) {
+
 	codec2_encode(codec2, bits, buf);
-	fwrite(bits, sizeof(char), nbyte, fout);
+
+	if (unpacked) {
+            /* unpack bits, MSB first */
+
+            bit = 7; byte = 0;
+            for(i=0; i<nbit; i++) {
+                unpacked_bits[i] = (bits[byte] >> bit) & 0x1;
+                bit--;
+                if (bit < 0) {
+                    bit = 7;
+                    byte++;
+                }
+            }
+            fwrite(unpacked_bits, sizeof(int), nbit, fout);
+        }
+        else
+            fwrite(bits, sizeof(char), nbyte, fout);
+
 	// if this is in a pipeline, we probably don't want the usual
         // buffering to occur
+
         if (fout == stdout) fflush(stdout);
         if (fin == stdin) fflush(stdin);
     }
@@ -112,6 +139,7 @@ int main(int argc, char *argv[])
 
     free(buf);
     free(bits);
+    free(unpacked_bits);
     fclose(fin);
     fclose(fout);
 
diff --git a/codec2-dev/src/codec2.c b/codec2-dev/src/codec2.c
index 7289f54f..acc62166 100644
--- a/codec2-dev/src/codec2.c
+++ b/codec2-dev/src/codec2.c
@@ -68,8 +68,8 @@ void codec2_encode_1300(struct CODEC2 *c2, unsigned char * bits, short speech[])
 void codec2_decode_1300(struct CODEC2 *c2, short speech[], const unsigned char * bits, float ber_est);
 void codec2_encode_1200(struct CODEC2 *c2, unsigned char * bits, short speech[]);
 void codec2_decode_1200(struct CODEC2 *c2, short speech[], const unsigned char * bits);
-void codec2_encode_650(struct CODEC2 *c2, unsigned char * bits, short speech[]);
-void codec2_decode_650(struct CODEC2 *c2, short speech[], const unsigned char * bits);
+void codec2_encode_700(struct CODEC2 *c2, unsigned char * bits, short speech[]);
+void codec2_decode_700(struct CODEC2 *c2, short speech[], const unsigned char * bits);
 static void ear_protection(float in_out[], int n);
 
 /*---------------------------------------------------------------------------*\
@@ -108,7 +108,7 @@ struct CODEC2 * CODEC2_WIN32SUPPORT codec2_create(int mode)
 	   (mode == CODEC2_MODE_1400) || 
 	   (mode == CODEC2_MODE_1300) || 
 	   (mode == CODEC2_MODE_1200) ||
-	   (mode == CODEC2_MODE_650) 
+	   (mode == CODEC2_MODE_700) 
 	   );
     c2->mode = mode;
     for(i=0; i<M; i++)
@@ -202,8 +202,8 @@ int CODEC2_WIN32SUPPORT codec2_bits_per_frame(struct CODEC2 *c2) {
 	return 52;
     if  (c2->mode == CODEC2_MODE_1200)
 	return 48;
-    if  (c2->mode == CODEC2_MODE_650)
-	return 26;
+    if  (c2->mode == CODEC2_MODE_700)
+	return 28;
 
     return 0; /* shouldn't get here */
 }
@@ -232,7 +232,7 @@ int CODEC2_WIN32SUPPORT codec2_samples_per_frame(struct CODEC2 *c2) {
 	return 320;
     if  (c2->mode == CODEC2_MODE_1200)
 	return 320;
-    if  (c2->mode == CODEC2_MODE_650)
+    if  (c2->mode == CODEC2_MODE_700)
 	return 320;
 
     return 0; /* shouldnt get here */
@@ -248,7 +248,7 @@ void CODEC2_WIN32SUPPORT codec2_encode(struct CODEC2 *c2, unsigned char *bits, s
 	   (c2->mode == CODEC2_MODE_1400) || 
 	   (c2->mode == CODEC2_MODE_1300) || 
 	   (c2->mode == CODEC2_MODE_1200) ||
-	   (c2->mode == CODEC2_MODE_650)
+	   (c2->mode == CODEC2_MODE_700)
 	   );
 
     if (c2->mode == CODEC2_MODE_3200)
@@ -263,8 +263,8 @@ void CODEC2_WIN32SUPPORT codec2_encode(struct CODEC2 *c2, unsigned char *bits, s
 	codec2_encode_1300(c2, bits, speech);
     if (c2->mode == CODEC2_MODE_1200)
 	codec2_encode_1200(c2, bits, speech);
-    if (c2->mode == CODEC2_MODE_650)
-	codec2_encode_650(c2, bits, speech);
+    if (c2->mode == CODEC2_MODE_700)
+	codec2_encode_700(c2, bits, speech);
 }
 
 void CODEC2_WIN32SUPPORT codec2_decode(struct CODEC2 *c2, short speech[], const unsigned char *bits)
@@ -282,7 +282,7 @@ void CODEC2_WIN32SUPPORT codec2_decode_ber(struct CODEC2 *c2, short speech[], co
 	   (c2->mode == CODEC2_MODE_1400) || 
 	   (c2->mode == CODEC2_MODE_1300) || 
 	   (c2->mode == CODEC2_MODE_1200) ||
-	   (c2->mode == CODEC2_MODE_650)
+	   (c2->mode == CODEC2_MODE_700)
 	   );
 
     if (c2->mode == CODEC2_MODE_3200)
@@ -297,8 +297,8 @@ void CODEC2_WIN32SUPPORT codec2_decode_ber(struct CODEC2 *c2, short speech[], co
  	codec2_decode_1300(c2, speech, bits, ber_est);
     if (c2->mode == CODEC2_MODE_1200)
  	codec2_decode_1200(c2, speech, bits);
-    if (c2->mode == CODEC2_MODE_650)
- 	codec2_decode_650(c2, speech, bits);
+    if (c2->mode == CODEC2_MODE_700)
+ 	codec2_decode_700(c2, speech, bits);
 }
 
 
@@ -1317,11 +1317,11 @@ void codec2_decode_1200(struct CODEC2 *c2, short speech[], const unsigned char *
 
 /*---------------------------------------------------------------------------*\
                                                        
-  FUNCTION....: codec2_encode_650	     
+  FUNCTION....: codec2_encode_700	     
   AUTHOR......: David Rowe			      
   DATE CREATED: April 2015
 
-  Encodes 320 speech samples (40ms of speech) into 26 bits.  
+  Encodes 320 speech samples (40ms of speech) into 28 bits.  
 
   The codec2 algorithm actually operates internally on 10ms (80
   sample) frames, so we run the encoding algorithm four times:
@@ -1329,7 +1329,7 @@ void codec2_decode_1200(struct CODEC2 *c2, short speech[], const unsigned char *
   frame 0: nothing
   frame 1: nothing
   frame 2: nothing
-  frame 3: voicing bit, scalar Wo and E, 17 bit LSP MEL scalar
+  frame 3: voicing bit, scalar Wo and E, 17 bit LSP MEL scalar, 2 spare
 
   The bit allocation is:
 
@@ -1339,11 +1339,12 @@ void codec2_decode_1200(struct CODEC2 *c2, short speech[], const unsigned char *
     Energy                              0          3         3
     log Wo                              0          5         5
     Voicing                             0          1         1
-    TOTAL                               0         26        26
+    spare                               0          2         2
+    TOTAL                               0         28        28
  
 \*---------------------------------------------------------------------------*/
 
-void codec2_encode_650(struct CODEC2 *c2, unsigned char * bits, short speech[])
+void codec2_encode_700(struct CODEC2 *c2, unsigned char * bits, short speech[])
 {
     MODEL   model;
     float   lsps[LPC_ORD_LOW];
@@ -1355,6 +1356,7 @@ void codec2_encode_650(struct CODEC2 *c2, unsigned char * bits, short speech[])
     unsigned int nbit = 0;
     float   bpf_out[4*N];
     short   bpf_speech[4*N];
+    int     spare = 0;
 
     assert(c2 != NULL);
 
@@ -1403,21 +1405,23 @@ void codec2_encode_650(struct CODEC2 *c2, unsigned char * bits, short speech[])
         pack(bits, &nbit, indexes[i], mel_bits(i));
     }
 
+    pack_natural_or_gray(bits, &nbit, spare, 2, c2->gray);
+
     assert(nbit == (unsigned)codec2_bits_per_frame(c2));
 }
 
 
 /*---------------------------------------------------------------------------*\
                                                        
-  FUNCTION....: codec2_decode_650	     
+  FUNCTION....: codec2_decode_700	     
   AUTHOR......: David Rowe			      
   DATE CREATED: April 2015
 
-  Decodes frames of 26 bits into 320 samples (40ms) of speech.
+  Decodes frames of 28 bits into 320 samples (40ms) of speech.
 
 \*---------------------------------------------------------------------------*/
 
-void codec2_decode_650(struct CODEC2 *c2, short speech[], const unsigned char * bits)
+void codec2_decode_700(struct CODEC2 *c2, short speech[], const unsigned char * bits)
 {
     MODEL   model[4];
     int     indexes[LPC_ORD_LOW];
diff --git a/codec2-dev/src/codec2.h b/codec2-dev/src/codec2.h
index 77c0f2dd..834d9a74 100644
--- a/codec2-dev/src/codec2.h
+++ b/codec2-dev/src/codec2.h
@@ -52,7 +52,7 @@
 #define CODEC2_MODE_1400 3
 #define CODEC2_MODE_1300 4
 #define CODEC2_MODE_1200 5
-#define CODEC2_MODE_650  6
+#define CODEC2_MODE_700  6
 
 struct CODEC2;
 
diff --git a/codec2-dev/src/cohpsk_ch.c b/codec2-dev/src/cohpsk_ch.c
index 78f36c1c..16c3791b 100644
--- a/codec2-dev/src/cohpsk_ch.c
+++ b/codec2-dev/src/cohpsk_ch.c
@@ -47,10 +47,12 @@
 #include "comp_prim.h"
 #include "../unittest/noise_samples.h"
 #include "ht_coeff.h"
+#include "ssbfilt_coeff.h"
 #include "codec2_fdmdv.h"
 
 #define BUF_N        160
 #define HF_DELAY_MS  2.0
+#define PAPR_TARGET  7.0
 
 /* This file gets generated using the function write_noise_file in tcohpsk.m.  You have to run
    tcohpsk first (any variant) to load the function into Octave, e.g.:
@@ -64,25 +66,27 @@
 int main(int argc, char *argv[])
 {
     FILE          *fin, *ffading, *fout;
-    float          EsNodB, foff_hz;
+    float          NodB, foff_hz;
     int            fading_en, nhfdelay;
 
     short          buf[BUF_N];
     float          htbuf[HT_N+BUF_N];
     COMP           ch_in[BUF_N];
     COMP           ch_fdm[BUF_N];
+    float          ssbfiltbuf[SSBFILT_N+BUF_N];
+    float          ssbfiltout[BUF_N];
                                            
     COMP           phase_ch;
     int            noise_r, noise_end;
-    float          EsNo, variance;
+    float          No, variance;
     COMP           scaled_noise;
     float          hf_gain;
     COMP          *ch_fdm_delay, aspread, aspread_2ms, delayed, direct;
-    float          tx_pwr, rx_pwr, noise_pwr;
-    int            frames, i, j, k, ret;
-    float          sam;
+    float          tx_pwr, tx_pwr_fade, noise_pwr;
+    int            frames, i, j, k, ret, clipped;
+    float          sam, peak, inclip, papr, CNo, snr3k, EbNo700;
 
-    if (argc == 6) {
+    if (argc == 7) {
         if (strcmp(argv[1], "-")  == 0) fin = stdin;
         else if ( (fin = fopen(argv[1],"rb")) == NULL ) {
             fprintf(stderr, "Error opening input modem raw file: %s: %s.\n",
@@ -97,28 +101,29 @@ int main(int argc, char *argv[])
             exit(1);
         }
 
-        EsNodB = atof(argv[3]);
+        NodB = atof(argv[3]);
         foff_hz = atof(argv[4]);
         fading_en = atoi(argv[5]);
+        inclip = atof(argv[6]);
     }
     else {
-        fprintf(stderr, "usage: %s InputRealModemRawFileFs7500Hz OutputRealModemRawFileFs7500Hz SNR3000Hz FoffHz FadingEn\n", argv[0]);
+        fprintf(stderr, "usage: %s InputRealModemRawFileFs7500Hz OutputRealModemRawFileFs7500Hz SNR3000Hz FoffHz FadingEn InputClip0to1\n", argv[0]);
         exit(1);
     }
-    fprintf(stderr, "EsNodB: %4.2f foff: %4.2f Hz fading: %d\n", EsNodB, foff_hz, fading_en);
+    fprintf(stderr, "NodB: %4.2f foff: %4.2f Hz fading: %d inclip: %4.2f\n", NodB, foff_hz, fading_en, inclip);
     
     phase_ch.real = 1.0; phase_ch.imag = 0.0; 
     noise_r = 0; 
     noise_end = sizeof(noise)/sizeof(COMP);
 
-    /*  each carrier has power = 2, total power 2Nc, total symbol rate
-        NcRs, noise BW B=Fs Es/No = (C/Rs)/(N/B), N = var =
-        2NcFs/NcRs(Es/No) = 2Fs/Rs(Es/No) */
+    /*  N = var = NoFs */
 
-    EsNo = pow(10.0, EsNodB/10.0);
-    variance = 2.0*COHPSK_FS/(COHPSK_RS*EsNo);
+    No = pow(10.0, NodB/10.0);
+    variance = COHPSK_FS*No;
     
-    tx_pwr = rx_pwr = noise_pwr = 0.0;
+    tx_pwr = tx_pwr_fade = noise_pwr = 0.0;
+    clipped = 0;
+    peak = 0.0;
 
     /* init HF fading model */
 
@@ -166,7 +171,19 @@ int main(int argc, char *argv[])
                freqencies.
             */
 
-            htbuf[j] = (float)buf[i]/FDMDV_SCALE;
+            sam = (float)buf[i];
+            //printf("sam: %f ", sam);
+            if (sam > inclip*32767.0)
+                sam = inclip*32767.0;
+            if (sam < -inclip*32767.0)
+                sam = -inclip*32767.0;
+            //printf("sam: %f\n", sam);
+            htbuf[j] = sam/FDMDV_SCALE;
+            
+            if (fabs(htbuf[j]) > peak) {
+                peak = fabs(htbuf[j]);
+            }
+            tx_pwr += pow(htbuf[j], 2.0);
 
             /* FIR filter with HT to get imag, just delay to get real */
 
@@ -185,19 +202,6 @@ int main(int argc, char *argv[])
         for(i=0; i<HT_N; i++)
            htbuf[i] = htbuf[i+BUF_N];
 
-        for(i=0; i<BUF_N; i++) {
-            //printf("%d %f %f\n", i, ch_in[i].real, ch_in[i].imag);
-            tx_pwr += pow(ch_in[i].real, 2.0) + pow(ch_in[i].imag, 2.0);
-        }
-
-        /* +3dB factor is beacuse we ouput a real signal, this has half
-           the power of the complex version but as the noise reflects
-           across from the -ve side the same noise power. */
-
-        for(i=0; i<BUF_N; i++) {
-            ch_in[i] = fcmult(sqrt(2.0), ch_in[i]);
-        }
-
 	/* --------------------------------------------------------*\
 	                          Channel
 	\*---------------------------------------------------------*/
@@ -231,10 +235,12 @@ int main(int argc, char *argv[])
             }
         }
 
-        /* we only output the real signal, which is half the power. */
+        /* Measure power after fading model to make sure aaverage pwr
+           is the same as AWGN channels. We only output the real
+           signal, which is half the power. */
 
         for(i=0; i<BUF_N; i++) {
-            rx_pwr += pow(ch_fdm[i].real, 2.0);
+            tx_pwr_fade += pow(ch_fdm[i].real, 2.0);
         }
 
         /* AWGN noise ------------------------------------------*/
@@ -242,21 +248,44 @@ int main(int argc, char *argv[])
         for(i=0; i<BUF_N; i++) {
             scaled_noise = fcmult(sqrt(variance), noise[noise_r]);
             ch_fdm[i] = cadd(ch_fdm[i], scaled_noise);
-            noise_pwr += pow(noise[noise_r].real, 2.0) + pow(noise[noise_r].imag, 2.0);
+            noise_pwr += pow(scaled_noise.real, 2.0);
             noise_r++;
             if (noise_r > noise_end) {
                 noise_r = 0;
                 //fprintf(stderr, "  [%d] noise wrap\n", f);            
+            }              
+        }
+
+        /* FIR filter to simulate (a rather flat) SSB filter.  Might
+           be useful to have an option for a filter with a few dB
+           ripple too, to screw up the modem. This is mainly so analog
+           SSB sounds realistic. */
+
+        for(i=0, j=SSBFILT_N; i<BUF_N; i++,j++) {
+            ssbfiltbuf[j] = ch_fdm[i].real;
+            ssbfiltout[i] = 0.0;
+            for(k=0; k<SSBFILT_N; k++) {
+                ssbfiltout[i] += ssbfiltbuf[j-k]*ssbfilt_coeff[k];
             }
-               
         }
 
+        /* update SSB filter memory */
+
+        for(i=0; i<SSBFILT_N; i++)
+           ssbfiltbuf[i] = ssbfiltbuf[i+BUF_N];
+
 	/* scale and save to disk as shorts */
 
 	for(i=0; i<BUF_N; i++) {
-            sam = FDMDV_SCALE * ch_fdm[i].real;
-            if (fabs(sam) > 32767.0)
-                fprintf(stderr,"clipping: %f\n", sam);
+            sam = FDMDV_SCALE * ssbfiltout[i];
+            if (sam > 32767.0) {
+                clipped++;
+                sam = 32767.0;
+            }
+            if (sam < -32767.0) {
+                clipped++;
+                sam = -32767.0;
+            }
 	    buf[i] = sam;
         }
 
@@ -271,13 +300,26 @@ int main(int argc, char *argv[])
 
     fclose(fin);
     fclose(fout);
-
-    fprintf(stderr, "tx var: %f noise var: %f rx var: %f\n", 
-           tx_pwr/(frames*BUF_N), 
-           noise_pwr/(frames*BUF_N),
-           rx_pwr/(frames*BUF_N) 
+    
+    fprintf(stderr, "peak pwr.....: %7.2f\nav input pwr.: %7.2f\nav pwr fading: %7.2f\nnoise pwr....: %7.2f\nclipping.....: %7.2f %%\n", 
+            peak*peak,
+            tx_pwr/(frames*BUF_N), 
+            tx_pwr_fade/(frames*BUF_N),
+            noise_pwr/(frames*BUF_N),
+            100.0*clipped/frames
            );
-
+    papr = 10*log10(peak*peak/(tx_pwr/(frames*BUF_N)));
+    CNo = 10*log10(tx_pwr/(noise_pwr/(COHPSK_FS/2))); // single sided spctrum magic IDFK!
+    snr3k = CNo - 10*log10(3000);
+    EbNo700 = CNo - 10*log10(700) - 10*log10(6.0/4.0); // divide by bit rate and pilot overhead
+    fprintf(stderr, "PAPR (dB)....: %7.2f (target %3.2f)\nC/No (dB)....: %7.2f\nSNR3k........: %7.2f\nEb/No(Rb=700): %7.2f\n",
+            papr,
+            PAPR_TARGET,
+            CNo,
+            snr3k,
+            EbNo700
+            );
+ 
     return 0;
 }
 
diff --git a/codec2-dev/src/cohpsk_demod.c b/codec2-dev/src/cohpsk_demod.c
index fdac1a94..507d06ac 100644
--- a/codec2-dev/src/cohpsk_demod.c
+++ b/codec2-dev/src/cohpsk_demod.c
@@ -34,20 +34,28 @@
 #include <errno.h>
 
 #include "codec2_cohpsk.h"
+#include "cohpsk_defs.h"
+#include "cohpsk_internal.h"
 #include "codec2_fdmdv.h"
+#include "octave.h"
+
+#define LOG_FRAMES 100
 
 int main(int argc, char *argv[])
 {
-    FILE          *fin, *fout;
+    FILE          *fin, *fout, *foct;
     struct COHPSK *cohpsk;
     int           rx_bits[COHPSK_BITS_PER_FRAME];
     COMP          rx_fdm[COHPSK_SAMPLES_PER_FRAME];
     short         rx_fdm_scaled[COHPSK_SAMPLES_PER_FRAME];
     int           frames, reliable_sync_bit, nin_frame;
-    int           i;
+    float        *rx_amp_log;
+    float        *rx_phi_log;
+    COMP         *rx_symb_log;
+    int           i, r, c, log_data_r, oct, logframes;
 
     if (argc < 3) {
-	printf("usage: %s InputModemRawFile OutputOneBitPerIntFile\n", argv[0]);
+	printf("usage: %s InputModemRawFile OutputOneBitPerIntFile [OctaveLogFile]\n", argv[0]);
 	exit(1);
     }
 
@@ -65,8 +73,30 @@ int main(int argc, char *argv[])
 	exit(1);
     }
 
+    foct = NULL;
+    oct = 0;
+    if (argc == 4) {
+        if ( (foct = fopen(argv[3],"wt")) == NULL ) {
+            fprintf(stderr, "Error opening output Octave file: %s: %s.\n",
+                    argv[3], strerror(errno));            
+	exit(1);
+        }
+        oct = 1;
+    }
+
     cohpsk = cohpsk_create();
 
+    if (oct) {
+        logframes = LOG_FRAMES;
+        rx_amp_log = (float *)malloc(sizeof(float)*logframes*NSYMROW*COHPSK_NC*ND);
+        assert(rx_amp_log != NULL);
+        rx_phi_log = (float *)malloc(sizeof(float)*logframes*NSYMROW*COHPSK_NC*ND);
+        assert(rx_phi_log != NULL);
+        rx_symb_log = (COMP *)malloc(sizeof(COMP)*logframes*NSYMROW*COHPSK_NC*ND);
+        assert(rx_symb_log != NULL);
+    }
+
+    log_data_r = 0;
     frames = 0;
 
     nin_frame = COHPSK_SAMPLES_PER_FRAME;
@@ -82,9 +112,23 @@ int main(int argc, char *argv[])
 
 	cohpsk_demod(cohpsk, rx_bits, &reliable_sync_bit, rx_fdm, &nin_frame);
 
- 	if (reliable_sync_bit)
+ 	if (reliable_sync_bit) {
             fwrite(rx_bits, sizeof(int), COHPSK_BITS_PER_FRAME, fout);
 
+            if (oct) {
+                for(r=0; r<NSYMROW; r++, log_data_r++) {
+                    for(c=0; c<COHPSK_NC*ND; c++) {
+                        rx_amp_log[log_data_r*COHPSK_NC*ND+c] = cohpsk->amp_[r][c]; 
+                        rx_phi_log[log_data_r*COHPSK_NC*ND+c] = cohpsk->phi_[r][c]; 
+                        rx_symb_log[log_data_r*COHPSK_NC*ND+c] = cohpsk->rx_symb[r][c]; 
+                    }
+                }
+                //printf("frames: %d log_data_r: %d\n", frames, log_data_r);
+                if (frames == logframes)
+                    oct = 0;
+            }
+        }
+
 	/* if this is in a pipeline, we probably don't want the usual
 	   buffering to occur */
 
@@ -94,7 +138,18 @@ int main(int argc, char *argv[])
 
     fclose(fin);
     fclose(fout);
+
+    /* optionally dump Octave files */
+
+    if (foct != NULL) {
+        octave_save_float(foct, "rx_amp_log_c", (float*)rx_amp_log, log_data_r, COHPSK_NC*ND, COHPSK_NC*ND);  
+        octave_save_float(foct, "rx_phi_log_c", (float*)rx_phi_log, log_data_r, COHPSK_NC*ND, COHPSK_NC*ND);  
+        octave_save_complex(foct, "rx_symb_log_c", (COMP*)rx_symb_log, log_data_r, COHPSK_NC*ND, COHPSK_NC*ND);  
+        fclose(foct);
+    }
+
     cohpsk_destroy(cohpsk);
 
+ 
     return 0;
 }
diff --git a/codec2-dev/src/ssbfilt_coeff.h b/codec2-dev/src/ssbfilt_coeff.h
new file mode 100644
index 00000000..457f03a8
--- /dev/null
+++ b/codec2-dev/src/ssbfilt_coeff.h
@@ -0,0 +1,107 @@
+/* 600 - 2600 Hz FIR filter coeffs */
+/* Generated by make_ssbfilt Octave script */
+
+#define SSBFILT_N 100
+
+float ssbfilt_coeff[]={
+  0.000065,
+  -0.000030,
+  0.000041,
+  0.000010,
+  -0.000128,
+  -0.000072,
+  -0.000007,
+  -0.000095,
+  -0.000063,
+  0.000016,
+  -0.000000,
+  0.000022,
+  -0.000115,
+  -0.000233,
+  -0.000023,
+  -0.000315,
+  -0.000725,
+  0.000073,
+  0.000380,
+  -0.000345,
+  0.000895,
+  0.002401,
+  0.001241,
+  0.001409,
+  0.003106,
+  0.001236,
+  -0.001117,
+  -0.001091,
+  -0.003184,
+  -0.005981,
+  -0.006904,
+  -0.007920,
+  -0.005588,
+  -0.002546,
+  -0.003476,
+  0.005155,
+  0.017465,
+  0.010772,
+  0.013033,
+  0.035082,
+  0.018466,
+  -0.010261,
+  0.016676,
+  0.004890,
+  -0.076807,
+  -0.055969,
+  -0.007360,
+  -0.155769,
+  -0.203150,
+  0.179458,
+  0.475523,
+  0.179458,
+  -0.203150,
+  -0.155769,
+  -0.007360,
+  -0.055969,
+  -0.076807,
+  0.004890,
+  0.016676,
+  -0.010261,
+  0.018466,
+  0.035082,
+  0.013033,
+  0.010772,
+  0.017465,
+  0.005155,
+  -0.003476,
+  -0.002546,
+  -0.005588,
+  -0.007920,
+  -0.006904,
+  -0.005981,
+  -0.003184,
+  -0.001091,
+  -0.001117,
+  0.001236,
+  0.003106,
+  0.001409,
+  0.001241,
+  0.002401,
+  0.000895,
+  -0.000345,
+  0.000380,
+  0.000073,
+  -0.000725,
+  -0.000315,
+  -0.000023,
+  -0.000233,
+  -0.000115,
+  0.000022,
+  -0.000000,
+  0.000016,
+  -0.000063,
+  -0.000095,
+  -0.000007,
+  -0.000072,
+  -0.000128,
+  0.000010,
+  0.000041,
+  -0.000030
+};
\ No newline at end of file