Merge pull request drowe67#285 from drowe67/dr-sinusoidal

Sinusoidal model scripts maintenance

Author: drowe67

octave/plamp.m

@@ -82,7 +82,7 @@ function plamp(samname, f, samname2)
     clf;
     s = [ Sn(2*f-1,:) Sn(2*f,:) ];
     plot(s,'b');
-    if (file_in_path(".",model2_name))
+    if (nargin == 3)
       s2 = [ Sn2(2*f-1,:) Sn2(2*f,:) ];
       hold on; plot(s2,'r'); hold off;
     end
@@ -97,7 +97,7 @@ function plamp(samname, f, samname2)
     hold on;
     if plot_sw; plot((0:255)*4000/256, Sw(f,:),";Sw;b"); end
 
-    if (file_in_path(".",model2_name))
+    if (nargin == 3)
       Wo2 = model2(f,1);
       L2 = model2(f,2);
       Am2 = model2(f,3:(L2+2));

src/CMakeLists.txt

@@ -33,10 +33,10 @@ set(CODEBOOKSD
 )
 
 # LSP VQ (1200 mode)
-set(CODEBOOKSJVM
-    ${D}/lspjvm1.txt
-    ${D}/lspjvm2.txt
-    ${D}/lspjvm3.txt
+set(CODEBOOKSJMV
+    ${D}/lspjmv1.txt
+    ${D}/lspjmv2.txt
+    ${D}/lspjmv3.txt
 )
 
 # Joint Wo-energy VQ (1200 mode)
@@ -102,11 +102,11 @@ add_custom_command(
     DEPENDS generate_codebook ${CODEBOOKSD}
 )
 
-# codebookjvm.c
+# codebookjmv.c
 add_custom_command(
-    OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/codebookjvm.c
-    COMMAND generate_codebook lsp_cbjvm ${CODEBOOKSJVM} > ${CMAKE_CURRENT_BINARY_DIR}/codebookjvm.c
-    DEPENDS generate_codebook ${CODEBOOKSJVM}
+    OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/codebookjmv.c
+    COMMAND generate_codebook lsp_cbjmv ${CODEBOOKSJMV} > ${CMAKE_CURRENT_BINARY_DIR}/codebookjmv.c
+    DEPENDS generate_codebook ${CODEBOOKSJMV}
 )
 
 # codebookge.c
@@ -177,7 +177,7 @@ set(CODEC2_SRCS
     pack.c
     codebook.c
     codebookd.c
-    codebookjvm.c
+    codebookjmv.c
     codebookge.c
     codebooknewamp1.c
     codebooknewamp1_energy.c

src/c2sim.c

@@ -50,6 +50,7 @@
 #include "bpfb.h"
 #include "newamp1.h"
 #include "lpcnet_freq.h"
+#include "sd.h"
 
 void synth_one_frame(int n_samp, codec2_fftr_cfg fftr_inv_cfg, short buf[], MODEL *model, float Sn_[], float Pn[], int prede, float *de_mem, float gain);
 void print_help(const struct option *long_options, int num_opts, char* argv[]);
@@ -70,7 +71,7 @@ int main(int argc, char *argv[])
 
     int lpc_model = 0, order = LPC_ORD;
     int lsp = 0, lspd = 0, lspvq = 0;
-    int lspjvm = 0;
+    int lspjmv = 0;
     int prede = 0;
     int   postfilt;
     int   hand_voicing = 0, hi = 0, simlpcpf = 0, modelin=0, modelout=0;
@@ -92,7 +93,7 @@ int main(int argc, char *argv[])
     FILE *fam = NULL, *fWo = NULL;
     FILE *faw = NULL;
     FILE *fhm = NULL;
-    FILE *fjvm = NULL;
+    FILE *fjmv = NULL;
     FILE *flspEWov = NULL;
     FILE *ften_ms_centre = NULL;
     FILE *fmodelout = NULL;
@@ -133,7 +134,7 @@ int main(int argc, char *argv[])
         { "lsp", no_argument, &lsp, 1 },
         { "lspd", no_argument, &lspd, 1 },
         { "lspvq", no_argument, &lspvq, 1 },
-        { "lspjvm", no_argument, &lspjvm, 1 },
+        { "lspjmv", no_argument, &lspjmv, 1 },
         { "phase0", no_argument, &phase0, 1 },
         { "postfilter", no_argument, &postfilt, 1 },
         { "hand_voicing", required_argument, &hand_voicing, 1 },
@@ -282,7 +283,7 @@ int main(int argc, char *argv[])
                     exit(1);
                 }
 	    } else if(strcmp(long_options[option_index].name, "dump_pitch_e") == 0) {
-	        if ((fjvm = fopen(optarg,"wt")) == NULL) {
+	        if ((fjmv = fopen(optarg,"wt")) == NULL) {
 	            fprintf(stderr, "Error opening pitch & energy dump file: %s: %s.\n",
 		        optarg, strerror(errno));
                     exit(1);
@@ -399,7 +400,7 @@ int main(int argc, char *argv[])
                } else if(strcmp(optarg,"1200") == 0) {
 	            lpc_model = 1;
 		    scalar_quant_Wo_e = 1;
-	            lspjvm = 1;
+	            lspjmv = 1;
 	            phase0 = 1;
 	            postfilt = 1;
 	            decimate = 4;
@@ -497,7 +498,8 @@ int main(int argc, char *argv[])
     COMP Aw[FFT_ENC];
     COMP H[MAX_AMP];
 
-
+    float sd_sum = 0.0; int sd_frames = 0;
+    
     for(i=0; i<m_pitch; i++) {
 	Sn[i] = 1.0;
 	Sn_pre[i] = 1.0;
@@ -677,7 +679,7 @@ int main(int argc, char *argv[])
             #endif
 
 	    if (dump_pitch_e)
-		fprintf(fjvm, "%f %f %d ", model.Wo, snr, model.voiced);
+		fprintf(fjmv, "%f %f %d ", model.Wo, snr, model.voiced);
 
             #ifdef DUMP
 	    dump_snr(snr);
@@ -701,6 +703,7 @@ int main(int argc, char *argv[])
 	\*------------------------------------------------------------*/
 
 	if (lpc_model) {
+            float ak_[LPC_ORD+1];
 
             e = speech_to_uq_lsps(lsps, ak, Sn, w, m_pitch, order);
             for(i=0; i<order; i++)
@@ -712,7 +715,7 @@ int main(int argc, char *argv[])
             #endif
 
 	    if (dump_pitch_e)
-		fprintf(fjvm, "%f\n", e);
+		fprintf(fjmv, "%f\n", e);
 
             #ifdef DUMP
             dump_lsp(lsps);
@@ -724,26 +727,33 @@ int main(int argc, char *argv[])
 		encode_lsps_scalar(lsp_indexes, lsps, LPC_ORD);
 		decode_lsps_scalar(lsps_, lsp_indexes, LPC_ORD);
 		bw_expand_lsps(lsps_, LPC_ORD, 50.0, 100.0);
-		lsp_to_lpc(lsps_, ak, LPC_ORD);
+		lsp_to_lpc(lsps_, ak_, LPC_ORD);
 	    }
 
 	    if (lspd) {
 		encode_lspds_scalar(lsp_indexes, lsps, LPC_ORD);
 		decode_lspds_scalar(lsps_, lsp_indexes, LPC_ORD);
-		lsp_to_lpc(lsps_, ak, LPC_ORD);
+		lsp_to_lpc(lsps_, ak_, LPC_ORD);
 	    }
 
-	    if (lspjvm) {
+	    if (lspjmv) {
 		/* Jean-Marc's multi-stage, split VQ */
-		lspjvm_quantise(lsps, lsps_, LPC_ORD);
+		lspjmv_quantise(lsps, lsps_, LPC_ORD);
 		{
 		    float lsps_bw[LPC_ORD];
 		    memcpy(lsps_bw, lsps_, sizeof(float)*order);
 		    bw_expand_lsps(lsps_bw, LPC_ORD, 50.0, 100.0);
-		    lsp_to_lpc(lsps_bw, ak, LPC_ORD);
+		    lsp_to_lpc(lsps_bw, ak_, LPC_ORD);
 		}
 	    }
 
+            if (lsp || lspd || lspjmv) {
+                sd_sum += spectral_dist(ak, ak_, LPC_ORD, fft_fwd_cfg, FFT_ENC);
+                sd_frames ++;
+            }
+
+            memcpy(ak, ak_, (LPC_ORD+1)*sizeof(float));
+
 	    if (scalar_quant_Wo_e) {
 		e = decode_energy(encode_energy(e, E_BITS), E_BITS);
                 model.Wo = decode_Wo(&c2const, encode_Wo(&c2const, model.Wo, WO_BITS), WO_BITS);
@@ -1068,7 +1078,9 @@ int main(int argc, char *argv[])
 	fclose(fout);
 
     if (lpc_model) {
-    	fprintf(stderr, "SNR av = %5.2f dB\n", sum_snr/frames);
+    	fprintf(stderr, "LPC->{Am} SNR av: %5.2f dB over %d frames\n", sum_snr/frames, frames);
+        if (lsp || lspd || lspjmv)
+            fprintf(stderr, "LSP quantiser SD: %5.2f dB*dB over %d frames\n", sd_sum/sd_frames, sd_frames);     
     }
     if (newamp1vq) {
     	fprintf(stderr, "var: %3.2f dB*dB\n", se/nse);
@@ -1087,7 +1099,7 @@ int main(int argc, char *argv[])
     if (fWo     != NULL) fclose(fWo);
     if (faw     != NULL) fclose(faw);
     if (fhm     != NULL) fclose(fhm);
-    if (fjvm    != NULL) fclose(fjvm);
+    if (fjmv    != NULL) fclose(fjmv);
     if (flspEWov != NULL) fclose(flspEWov);
     if (fphasenn != NULL) fclose(fphasenn);
     if (frateK != NULL) fclose(frateK);

src/ch.c

@@ -4,7 +4,7 @@
   AUTHOR......: David Rowe
   DATE CREATED: May 2015
 
-  Channel impairment program for testing command line versions of modems.
+  Channel simulation program for testing command line versions of modems.
 
 \*---------------------------------------------------------------------------*/
 
@@ -92,15 +92,17 @@ int main(int argc, char *argv[])
 
     if (argc < 3) {
     helpmsg:
-        fprintf(stderr, "usage: %s InputRealModemRawFile OutputRealModemRawFile [Options]\n"
+        fprintf(stderr, "Command line channel simulation tool.\n"
+                        "\n"
+                        "usage: %s InputRealModemRawFile OutputRealModemRawFile [Options]\n"
                         "\n"
                         "  real int16 input -> Gain -> Hilbert Transform -> clipper -> freq shift ->\n"
                         "  Multipath -> AWGN noise -> SSB filter -> real int16 output\n"
                         "\n"
                         "[--clip int16]         Hilbert clipper (clip complex signal magnitude, default 32767)\n"
-                        "[--complexout]         Optional int16 IQ complex output (default real)\n"
+                        "[--complexout]         Optional int16 IQ complex output (default real int16)\n"
                         "[--ctest]              Check PAPR is around 0dB, used to support ctests\n"
-                        "[--freqq FoffHz]       Frequency offset (default 0Hz)\n"
+                        "[--freq FoffHz]        Frequency offset (default 0Hz)\n"
                         "[--fading_dir Path]    path to multipath fading files (default 'unittest')\n"
                         "[--Fs SampleRateHz]    Sample rate of simulation (default 8000 Hz)\n"
                         "[--gain G]             Linear gain (default 1.0)\n"
@@ -301,11 +303,11 @@ int main(int argc, char *argv[])
         frames++;
 
         /* Hilbert Transform to produce complex signal so we can do
-           single sided freq shifts, HF channel modemsl, and analog compression.
-           Allows us to use real signal I/O which is handy.
+           single sided freq shifts, HF channel models, and analog
+           compression.  Allows us to use real signal I/O.
 
-           As the real and imag filters both have unity gain, ch_in[] has twice
-           the power of the real input signal buf[].
+           As the real and imag filters both have unity gain, ch_in[]
+           has twice the power of the real input signal buf[].
         */
 
         for(i=0, j=HT_N; i<BUF_N; i++,j++) {
@@ -363,7 +365,7 @@ int main(int argc, char *argv[])
                 ch_fdm_delay[i] = ch_fdm[j];
 
             /* combine direct and delayed paths, both multiplied by
-               "spreading" (doppler) functions */
+               "spreading" (Doppler) functions */
 
             for(i=0; i<BUF_N; i++) {
                 ret = fread(&aspread, sizeof(COMP), 1, ffading);
@@ -465,7 +467,7 @@ int main(int argc, char *argv[])
     if (ffading != NULL) fclose(ffading);
     if (ch_fdm_delay != NULL) FREE(ch_fdm_delay);
     if (ctest) {
-        /* special ctest modes, check CPAPR is around 0dB */
+        /* special ctest mode: check CPAPR is around 0dB */
         if (fabs(papr) < 0.7) return 0; else return 1;
     }
     else return 0;