anti Git rant and some mods to c2sim to support playing with mel-lsps for a new impro...
authordrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Wed, 5 Aug 2015 06:05:46 +0000 (06:05 +0000)
committerdrowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Wed, 5 Aug 2015 06:05:46 +0000 (06:05 +0000)
git-svn-id: https://svn.code.sf.net/p/freetel/code@2257 01035d8c-6547-0410-b346-abe4f91aad63

codec2-dev/README
codec2-dev/octave/cohpsk_frame_design.ods
codec2-dev/octave/melvq.m [new file with mode: 0644]
codec2-dev/octave/sd.m
codec2-dev/src/c2sim.c
codec2-dev/src/codec2.c
codec2-dev/src/dump.c
codec2-dev/src/dump.h

index a6e8b6055ee00c163c31812aeaad35c62fef6269..4cd93a59303883f5db79d603ed62ce8e32bc3517 100644 (file)
@@ -1,7 +1,7 @@
 Codec 2 README
 --------------
 
-Codec 2 is an open source (LGPL licensed) speech codec for 2400 bit/s
+Codec 2 is an open source (LGPL licensed) speech codec for 3200 bit/s
 and below.  For more information please see:
 
     http://rowetel.com/codec2.html
@@ -10,6 +10,25 @@ Also included is a FDMDV modem (README_fdmdv.txt), a coherent PSK
 modem and an API for embedding FreeDV in other programs (see example
 below).  For more information on building Codec 2 see README.cmake
 
+SVN Repository
+--------------
+
+Check out the latest (development branch) code using:
+
+  $ svn co https://svn.code.sf.net/p/freetel/code/codec2-dev codec2-dev
+
+There are 3rd party GIT mirrors of Codec2 and FreeDV. Use Git at your
+own risk.
+
+  GIT IS NOT SUPPORTED!!!
+
+All patches, support questions etc, need to be against the SVN
+repository above.
+
+Please do not email me (David), or the codec2-dev mailing list
+suggesting we change to Git.  I get these emails every week.  Really,
+I understand the arguments, but am content with SVN for now.
+
 Quickstart
 ----------
 
index f1b2c9187fddea7be3335a34c6547ef0fc8d0c6e..91d16ad491bb21c95ff24e5ce09e84e9e8299a60 100644 (file)
Binary files a/codec2-dev/octave/cohpsk_frame_design.ods and b/codec2-dev/octave/cohpsk_frame_design.ods differ
diff --git a/codec2-dev/octave/melvq.m b/codec2-dev/octave/melvq.m
new file mode 100644 (file)
index 0000000..aab2203
--- /dev/null
@@ -0,0 +1,130 @@
+% melvq.m
+% David Rowe Aug 2015
+%
+% Experimenting with VQ design for mel LSPs
+
+% todo:
+% [ ] Sorting inside search what happens if we get a order issue, fix and calc SD
+% [ ] Search to avoid a worst case error rather than average?  This could be included
+%     in training.
+% [ ] nested/staged search
+
+1;
+
+% train up multi-stage VQ
+
+function vq = trainvq(training_data, Nvec, stages)
+
+  vq = [];
+  for i=1:stages
+    [idx centers] = kmeans(training_data, Nvec);
+    quant_error = centers(idx,:) - training_data;
+    printf("mse stage %d: %f\n", i, mean(std(quant_error)));
+    training_data = quant_error;
+    vq(:,:,i) = centers;
+  end
+
+end
+
+
+function [mse_list index_list] = search_vq(vq, target, m)
+
+  [Nvec order] = size(vq);
+
+  mse = zeros(1, Nvec);
+  % find mse for each vector
+
+  for i=1:Nvec
+     mse(i) = sum((target - vq(i,:)) .^2);
+  end
+
+  % sort and keep top m matches
+
+  [mse_list index_list ] = sort(mse);
+
+  mse_list = mse_list(1:m);
+  index_list = index_list(1:m);
+
+endfunction
+
+
+% Search multi-stage VQ, retaining m best candidates at each stage
+
+function [res ind] = mbest(vqset, input_vecs, m)
+
+  [Nvec order stages] = size(vqset);
+  [Ninput tmp] = size(input_vecs);
+
+  res = [];   % residual error after VQ
+  ind = [];   % index of vqs
+
+  for i=1:Ninput
+  
+    % first stage, find mbest candidates
+
+    [mse_list index_list] = search_vq(vqset(:,:,1), input_vecs(i,:), m);
+    cand_list = [mse_list' index_list'];
+    cand_list = sortrows(cand_list,1);
+
+    % subsequent stages ...........
+
+    for s=2:stages
+
+      % compute m targets for next stage, and update path
+
+      prev_indexes = zeros(m,s-1);
+      for t=1:m
+        target(t,:) = input_vecs(i,:);
+        for v=1:s-1
+          target(t,:) -= vqset(cand_list(t,v+1),:,v);
+        end
+        prev_indexes(t,:) = cand_list(t,2:s);
+      end
+      
+      % search stage s using m targets from stage s-1
+      % with m targets, we do m searches which return the m best possibilities
+      % so we get a matrix with one row per candidate, m*m rows total
+      % prev_indexes provides us with the path through the VQs for each candidate row
+
+      avq = vqset(:,:,s);
+      cand_list = [];
+      for t=1:m
+        [mse_list index_list] = search_vq(avq, target(t,:), m);
+        x = ones(m,1)*prev_indexes(t,:);
+        cand_row = [mse_list' x index_list'];
+        cand_list = [cand_list; cand_row];
+      end
+
+      % sort into m best rows
+     
+      cand_list = sortrows(cand_list,1);
+      cand_list = cand_list(1:m,:);
+
+    end
+
+    % final residual
+    target(1,:) = input_vecs(i,:);
+    for v=1:s
+      target(1,:) -= vqset(cand_list(1,v+1),:,v);
+    end
+    res = [res; target(1,:)];
+    ind = [ind; cand_list(1,2:1+stages)];
+  end
+
+endfunction
+
+more off;
+load "../build_linux/src/all_mel.txt"
+load vq;
+[res1 ind] = mbest(vq, all_mel(1,:),3);
+mean(std(res1))
+
+% save text file of "vq quantised mels"
+% load back into c2sim at run time
+% train on continuous mels
+% sorting/stability
+% see how it sounds
+% Goal is to get VQ sounding OK, similar to UQ at 20 or 40ms dec,
+% sig better than current 700
+% analysis of data, outliers, different training
index 7995a9c2ce8974b2ec847aa1b98446885f45e8a0..d300524092d9f3787eac4e925cacbd9f95857074 100644 (file)
@@ -52,7 +52,7 @@ function sd(raw_filename, dump_file_prefix, f)
   axis([1 length(s) -20E3 20E3])
   subplot(212)
   [a b c] = plotyy(1:frames, sd, 1:frames, e);
-  axis(a, [1 frames 0 10])
+  %axis(a, [1 frames 0 10])
 
   lsp1_filename = sprintf("%s_lsp.txt", dump_file_prefix);
   lsp2_filename = sprintf("%s_lsp_.txt", dump_file_prefix);
@@ -77,9 +77,25 @@ function sd(raw_filename, dump_file_prefix, f)
   subplot(211)
   ind = find(e > 0);
   hist(sd(ind))
+  title('Histogram of SD');
   subplot(212)
   plot((1:Ndft)*8000/Ndft, spec_err)
   axis([300 3000 0 max(spec_err)])
+  title('Average error across spectrum')
+
+  mel_indexes_filename = sprintf("%s_mel_indexes.txt", dump_file_prefix);
+  if file_in_path(".", mel_indexes_filename)
+    mel_indexes = load(mel_indexes_filename);
+    figure(6)
+    bins = [15, 7, 15, 7, 7, 7];
+    ind = find(e > 5); % ignore silence frames
+    for i=1:6
+      subplot(3,2,i)
+      hist(mel_indexes(ind,i),0:bins(i))
+      ylab = sprintf("index %d", i);
+      ylabel(ylab);
+    end
+  end
 
   % now enter single step mode so we can analyse each frame
   k = ' ';
@@ -101,10 +117,10 @@ function sd(raw_filename, dump_file_prefix, f)
     figure(3);
     clf;
 
-    plot((1:Ndft/2)*4000/(Ndft/2), A1(fr,1:(Ndft/2)),";A1;r");
+    plot((1:Ndft/2)*4000/(Ndft/2), A1(fr,1:(Ndft/2)),";A enc;r");
     axis([1 4000 -20 40]);
     hold on;
-    plot((1:Ndft/2)*4000/(Ndft/2), A2(fr,1:(Ndft/2)),";A2;");
+    plot((1:Ndft/2)*4000/(Ndft/2), A2(fr,1:(Ndft/2)),";A dec;");
     if file_in_path(".",weights_filename)
       plot(lsp1(fr,:)*4000/pi, weights(fr,:),";weights;g+");
     end
@@ -118,9 +134,10 @@ function sd(raw_filename, dump_file_prefix, f)
     endfor
     printf("\n");
     
-    plot(0,0,';lsp1;r');
-    plot(0,0,';lsp2;b');
-    sd_str = sprintf(";sd %3.2f dB*dB;", sd(fr));
+    plot(0,0,';lsp enc;r');
+    plot(0,0,';lsp dec;b');
+    plot(0,0,';mel dec;g');
+    sd_str = sprintf(";sd %3.2f dB*dB;", sd(f));
     plot(0,0,sd_str);
    
     hold off;
index f0ed0efb9e25a41035175d71a247ddcfd5301787..89afbd04bdd5f3363d229bccfb5faedd995ab8d0 100644 (file)
@@ -235,8 +235,8 @@ int main(int argc, char *argv[])
             } else if(strcmp(long_options[option_index].name, "dec") == 0) {
                
                 decimate = atoi(optarg);
-               if ((decimate != 2) && (decimate != 4)) {
-                   fprintf(stderr, "Error in --dec, must be 2 or 4\n");
+               if ((decimate != 2) && (decimate != 3) && (decimate != 4)) {
+                   fprintf(stderr, "Error in --dec, must be 2, 3, or 4\n");
                    exit(1);
                }
 
@@ -624,20 +624,42 @@ int main(int argc, char *argv[])
                    mel[i] = floor(2595.0*log10(1.0 + f/700.0) + 0.5);
                }
 
+                /* round to the nearest x hz */
+                #define MEL_ROUND 50
+               for(i=0; i<order; i++) {
+                    float x = floor(mel[i]/MEL_ROUND+0.5)*MEL_ROUND;
+                    mel[i] = x;
+                    //printf("mel[%d] = %f x: %f\n", i, mel[i], x);
+                }
+
                for(i=1; i<order; i++) {
-                   if (mel[i] == mel[i-1])
-                       mel[i]++;
+                   if (mel[i] == mel[i-1]) {
+                       mel[i]+=MEL_ROUND/2;
+                       mel[i-1]-=MEL_ROUND/2;
+                        i = 1;
+                    }
                }
 
                encode_mels_scalar(mel_indexes, mel, 6);
+                #ifdef DUMP
+                dump_mel_indexes(mel_indexes, 6);
+                #endif
                decode_mels_scalar(mel, mel_indexes, 6);
                 
+               for(i=1; i<order; i++) {
+                   if (mel[i] == mel[i-1]) {
+                       mel[i]+=MEL_ROUND/2;
+                       mel[i-1]-=MEL_ROUND/2;
+                        i = 1;
+                    }
+               }
+
                 #ifdef DUMP
                 dump_mel(mel, order);
                 #endif
 
                for(i=0; i<LPC_ORD; i++) {
-                   f_ = 700.0*( pow(10.0, (float)mel[i]/2595.0) - 1.0);
+                   f_ = 700.0*( pow(10.0, mel[i]/2595.0) - 1.0);
                    lsps_[i] = f_*(PI/4000.0);
                }
  
index 8482fb87e5ae304fa0698622d492052bdc529cd1..c2eff53397885835f7e5c06e9af44da805c58574 100644 (file)
@@ -1372,7 +1372,7 @@ void codec2_encode_700(struct CODEC2 *c2, unsigned char * bits, short speech[])
         c2->bpf_buf[BPF_N+i] = speech[i];
     inverse_filter(&c2->bpf_buf[BPF_N], bpf, 4*N, bpf_out, BPF_N);
     for(i=0; i<4*N; i++)
-        bpf_speech[i] = bpf_out[i];
+        bpf_speech[i] = speech[i];
 
     /* frame 1 --------------------------------------------------------*/
 
index 826007270c5a7bae801255ce9a6371d6f73c128e..f172424d5badbddd3fc205c13ce3c686fcfe68aa 100644 (file)
@@ -56,6 +56,7 @@ static FILE *flsp = NULL;
 static FILE *fweights = NULL;
 static FILE *flsp_ = NULL;
 static FILE *fmel = NULL;
+static FILE *fmel_indexes = NULL;
 static FILE *fphase = NULL;
 static FILE *fphase_ = NULL;
 static FILE *ffw = NULL;
@@ -107,6 +108,8 @@ void dump_off(){
        fclose(flsp_);
     if (fmel != NULL)
        fclose(fmel);
+    if (fmel_indexes != NULL)
+       fclose(fmel_indexes);
     if (fphase != NULL)
        fclose(fphase);
     if (fphase_ != NULL)
@@ -494,6 +497,23 @@ void dump_mel(float mel[], int order) {
     fprintf(fmel,"\n");    
 }
 
+void dump_mel_indexes(int mel_indexes[], int order) {
+    int i;
+    char s[MAX_STR];
+
+    if (!dumpon) return;
+
+    if (fmel_indexes == NULL) {
+       sprintf(s,"%s_mel_indexes.txt", prefix);
+       fmel_indexes = fopen(s, "wt");
+       assert(fmel_indexes != NULL);
+    }
+
+    for(i=0; i<order; i++)
+       fprintf(fmel_indexes,"%d\t",mel_indexes[i]);
+    fprintf(fmel_indexes,"\n");    
+}
+
 void dump_ak(float ak[], int order) {
     int i;
     char s[MAX_STR];
index d816697c239efe2cc2053f22578e40240f81a2f7..3cf5868899da6b1270b51e80603905ba1965230c 100644 (file)
@@ -51,6 +51,7 @@ void dump_lsp(float lsp[]);
 void dump_weights(float w[], int ndim);
 void dump_lsp_(float lsp_[]);
 void dump_mel(float mel[], int order);
+void dump_mel_indexes(int mel_indexes[], int order);
 void dump_ak(float ak[], int order);
 void dump_ak_(float ak[], int order);
 void dump_E(float E);