plamp.m displays SNR on plot, moved back to SNR measure as per thesis

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

diff --git a/codec2/README.txt b/codec2/README.txt
index 27b19bec8bb20c2b6a8af49389a8c0b2ee590fd8..47b742611736b9f0614f696d88c3e5ff8ea41957 100644 (file)
--- a/codec2/README.txt
+++ b/codec2/README.txt
@@ -36,6 +36,32 @@ Plan
 [ ] Milestone 2 - Spectral amplitudes quantised 
 [ ] Milestone 3 - Prototype 2400 bit/s codec
 
+LPC Modelling
+-------------
+
+$ ./sinedec ../raw/hts1a.raw hts1a.mdl --lpc 10 - hts1a_lpc10.raw
+
+Disucss why LPC modelling works so well when Am recovered via RMS method
+(Section 5.1 of thesis).  Equal area model of LPC spectra versus harmonic?
+Seems to work remarkably well, especially compared to sampling.  SNRs up to
+30dB on female frames.
+
+Octave Scripts
+--------------
+
+pl.m    - plot a segment from a raw file
+
+pl2.m   - plot the same segments from two different files to compare
+
+plamp.m - menu based GUI interface to "dump" files, move back and forward
+          through file exaimining time and frequency domain parameters, lpc 
+          model etc
+  
+          $ ./sinedec ../raw/hts1a.raw hts1a.mdl --lpc 10 --dump hts1a
+          $ cd ../octave
+          $ octave
+          octave:1> plamp("../src/hts1",25)
+
 Directories
 -----------
 

diff --git a/codec2/octave/plamp.m b/codec2/octave/plamp.m
index 8ae51c6cedff13b7a36417fc04e749370b9f5a60..daaed16b18e23cd85eff4b6691b9484035a29e7d 100644 (file)
--- a/codec2/octave/plamp.m
+++ b/codec2/octave/plamp.m
@@ -1,11 +1,15 @@
 % Copyright David Rowe 2009
 % This program is distributed under the terms of the GNU General Public License 
 % Version 2
+%
+% Plot ampltiude modelling information from dump files.
 
 function plamp(samname, f)
   
   sn_name = strcat(samname,"_sn.txt");
   Sn = load(sn_name);
+  sw_name = strcat(samname,"_sw.txt");
+  Sw = load(sw_name);
   model_name = strcat(samname,"_model.txt");
   model = load(model_name);
   modelq_name = strcat(samname,"_qmodel.txt");
@@ -13,22 +17,42 @@ function plamp(samname, f)
   pw_name = strcat(samname,"_pw.txt");
   Pw = load(pw_name);
 
-  figure(1);
-  clg;
-  s = [ Sn(2*f+1,:) Sn(2*f+2,:) ];
-  plot(s);
- 
-  figure(2);
-  Wo = model(f,1);
-  L = model(f,2);
-  Am = model(f,3:(L+2));
-  plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;");
-  hold on;
-  Amq = modelq(f,3:(L+2));
-  plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;" );
-  plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;");
-  plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am),";Am_err;" );
-  hold off;
+  do 
+    figure(1);
+    clg;
+    s = [ Sn(2*f+1,:) Sn(2*f+2,:) ];
+    plot(s);
+    axis([1 length(s) -20000 20000]);
+
+    figure(2);
+    Wo = model(f,1);
+    L = model(f,2);
+    Am = model(f,3:(L+2));
+    plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;");
+    axis([1 4000 -10 80]);
+    hold on;
+    plot((0:255)*4000/256, Sw(f,:),";Sw;");
+    Amq = modelq(f,3:(L+2));
+    plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;" );
+    plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;");
+    signal = Am * Am';
+    noise = (Am-Amq) * (Am-Amq)'; 
+    snr = 10*log10(signal/noise);
+    Am_err_label = sprintf(";Am_err SNR %4.2f dB;",snr);
+    %plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am),";Am_err;" );
+    plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label);
+    hold off;
+
+    printf("\rframe: %d  menu: n-next  b-back  q-quit ", f);
+    fflush(stdout);
+    k = kbhit();
+    if (k == 'n')
+      f = f + 1;
+    endif
+    if (k == 'b')
+      f = f - 1;
+    endif
+  until (k == 'q')
+  printf("\n");
 
-  Am_err = mean(abs(20*log10(Amq) - 20*log10(Am)))
 endfunction

diff --git a/codec2/src/Makefile b/codec2/src/Makefile
index 164844b1e65593657efd33949318538def87f02f..b889ad04765050842cd9dbf762ccb32ff0df1c42 100644 (file)
--- a/codec2/src/Makefile
+++ b/codec2/src/Makefile
@@ -3,7 +3,7 @@ CFLAGS=-g -Wall
 
 SINENC_OBJ  = sinenc.o globals.o initenc.o four1.o refine.o spec.o
 SINEDEC_OBJ = sinedec.o globals.o initenc.o initdec.o four1.o synth.o \
-              quantise.o lpc.o dump.o
+              quantise.o lpc.o dump.o refine.o 
 
 all: sinenc sinedec
 

diff --git a/codec2/src/dump.c b/codec2/src/dump.c
index 0deb15cae5a07f8584826ef257a381c1bc4d04ec..5a8246765a7e78ef5fea2a0b93894d6deba39ac2 100644 (file)
--- a/codec2/src/dump.c
+++ b/codec2/src/dump.c
@@ -32,6 +32,7 @@
 static int dumpon = 0;
 
 static FILE *fsn;
+static FILE *fsw;
 static FILE *fmodel;
 static FILE *fqmodel;
 static FILE *fpw;
@@ -45,6 +46,10 @@ void dump_on(char prefix[]) {
     fsn = fopen(s, "wt");
     assert(fsn != NULL);
 
+    sprintf(s,"%s_sw.txt", prefix);
+    fsw = fopen(s, "wt");
+    assert(fsw != NULL);
+
     sprintf(s,"%s_model.txt", prefix);
     fmodel = fopen(s, "wt");
     assert(fmodel != NULL);
@@ -81,6 +86,17 @@ void dump_Sn(float Sn[]) {
     fprintf(fsn,"\n");    
 }
 
+void dump_Sw(COMP Sw[]) {
+    int i;
+
+    if (!dumpon) return;
+
+    for(i=0; i<FFT_ENC/2; i++)
+       fprintf(fsw,"%f\t",
+               10.0*log10(Sw[i].real*Sw[i].real + Sw[i].imag*Sw[i].imag));
+    fprintf(fsw,"\n");    
+}
+
 void dump_model(MODEL *model) {
     int l;
 

diff --git a/codec2/src/dump.h b/codec2/src/dump.h
index 61e41b2b1d7a6688216c2c92b2aa658d2a827eb4..8e4ccc61506bc962ebf1f49cf9245ea688de1fa1 100644 (file)
--- a/codec2/src/dump.h
+++ b/codec2/src/dump.h
@@ -32,6 +32,7 @@
 void dump_on(char filename_prefix[]);
 void dump_off();
 void dump_Sn(float Sn[]);
+void dump_Sw(COMP Sw[]);
 void dump_model(MODEL *m);
 void dump_quantised_model(MODEL *m);
 void dump_Pw(COMP Pw[]);

diff --git a/codec2/src/quantise.c b/codec2/src/quantise.c
index e995ba2c8a3767ab9ee1df6e9fac0ae4b81d9769..2a7e380feaba6ee49357e03a941e31a0a11451fa 100644 (file)
--- a/codec2/src/quantise.c
+++ b/codec2/src/quantise.c
@@ -54,7 +54,7 @@ float lpc_model_amplitudes(
   float ak[MAX_ORDER+1];
   float E;
   int   i;
-  float sd;                    /* spectral distortion for this frame */
+  float snr;   
 
   for(i=0; i<AW_ENC; i++)
       Wn[i] = Sn[i]*w[i];
@@ -64,9 +64,9 @@ float lpc_model_amplitudes(
   for(i=0; i<=order; i++)
       E += ak[i]*R[i];
   
-  aks_to_M2(ak,order,model,E,&sd);   /* {ak} -> {Am} LPC decode */
+  aks_to_M2(ak,order,model,E,&snr);   /* {ak} -> {Am} LPC decode */
 
-  return sd;
+  return snr;
 }
 
 /*---------------------------------------------------------------------------*\
@@ -84,7 +84,7 @@ void aks_to_M2(
   int   order,
   MODEL *model,        /* sinusoidal model parameters for this frame */
   float E,     /* energy term */
-  float *sd    /* spectral distortion for this frame in dB */
+  float *snr   /* signal to noise ratio for this frame in dB */
 )
 {
   COMP Pw[FFT_DEC];    /* power spectrum */
@@ -93,7 +93,7 @@ void aks_to_M2(
   float r;             /* no. rads/bin */
   float Em;            /* energy in band */
   float Am;            /* spectral amplitude sample */
-  float noise;
+  float signal, noise;
 
   r = TWO_PI/(FFT_DEC);
 
@@ -116,7 +116,7 @@ void aks_to_M2(
 
   /* Determine magnitudes by linear interpolation of P(w) -------------------*/
 
-  noise = 0.0;
+  signal = noise = 0.0;
   for(m=1; m<=model->L; m++) {
     am = floor((m - 0.5)*model->Wo/r + 0.5);
     bm = floor((m + 0.5)*model->Wo/r + 0.5);
@@ -126,9 +126,9 @@ void aks_to_M2(
       Em += Pw[i].real;
     Am = sqrt(Em);
 
-    noise += pow(log10(Am/model->A[m]),2.0);
+    signal += pow(model->A[m],2.0);
+    noise  += pow(model->A[m] - Am,2.0);
     model->A[m] = Am;
   }
-  *sd = 20.0*sqrt(noise/model->L);
-
+  *snr = 10.0*log10(signal/noise);
 }

diff --git a/codec2/src/quantise.h b/codec2/src/quantise.h
index 00583ab91f24e6661ff27d5c3af3204c609d8d76..5b1a15f0176fcd120ab0a418055909657ec0309b 100644 (file)
--- a/codec2/src/quantise.h
+++ b/codec2/src/quantise.h
@@ -30,6 +30,6 @@
 #include "sine.h"
 
 float lpc_model_amplitudes(float Sn[], MODEL *model, int order, int lsp);
-void aks_to_M2(float ak[], int order, MODEL *model, float E, float *sd);
+void aks_to_M2(float ak[], int order, MODEL *model, float E, float *snr);
 
 #endif

diff --git a/codec2/src/sinedec.c b/codec2/src/sinedec.c
index ea0de97d468969d1b19e8b6af8827c0bcbe39910..bd01343681ac87ab053a21ee5531b3c48e525749 100644 (file)
--- a/codec2/src/sinedec.c
+++ b/codec2/src/sinedec.c
@@ -72,8 +72,8 @@ int main(int argc, char *argv[])
 
   char  out_file[MAX_STR];
   int   arg;
-  float sd;
-  float sum_sd;
+  float snr;
+  float sum_snr;
 
   int lpc_model, order;
   int dump;
@@ -149,7 +149,7 @@ int main(int argc, char *argv[])
   /* Main loop ------------------------------------------------------------*/
 
   frames = 0;
-  sum_sd = 0;
+  sum_snr = 0;
   while(fread(&model,sizeof(model),1,fmodel)) {
     frames++;
 
@@ -161,14 +161,15 @@ int main(int argc, char *argv[])
     for(i=0; i<N; i++)
       Sn[i+N+AW_ENC/2] = buf[i];
     dump_Sn(Sn);
- 
+    dft_speech(); dump_Sw(Sw);   
+
     dump_model(&model);
 
     /* optional LPC model amplitudes */
 
     if (lpc_model) {
-       sd = lpc_model_amplitudes(Sn, &model, order, 0);
-       sum_sd += sd;
+       snr = lpc_model_amplitudes(Sn, &model, order, 0);
+       sum_snr += snr;
     }
 
     dump_quantised_model(&model);
@@ -197,7 +198,7 @@ int main(int argc, char *argv[])
     fclose(fout);
 
   if (lpc_model)
-      printf("sd av = %5.2f dB\n", sum_sd/frames);
+      printf("SNR av = %5.2f dB\n", sum_snr/frames);
 
   if (dump)
       dump_off();