new phase model and script to plot LSP PDFs

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

diff --git a/codec2/octave/lsp_pdf.m b/codec2/octave/lsp_pdf.m
new file mode 100644 (file)
index 0000000..a5a1b9a
--- /dev/null
+++ b/codec2/octave/lsp_pdf.m
@@ -0,0 +1,50 @@
+% lsp_pdf.m
+% David Rowe 2 Oct 2009
+% Plots histograms (PDF estimates) of LSP training data
+
+function lsp_pdf(lsp)
+  [r,c] = size(lsp);
+
+  % LSPs
+
+  figure(1);
+  clg;
+  [x,y] = hist(lsp(:,1),100);
+  plot(y,x,";1;");
+  hold on;
+  for i=2:c
+    [x,y] = hist(lsp(:,i),100);
+    legend = sprintf(";%d;",i);
+    plot(y,x,legend);
+  endfor
+  hold off;
+  grid;
+
+  % LSP differences
+
+  figure(2);
+  clg;
+  subplot(211)
+  [x,y] = hist(lsp(:,1),100);
+  plot(y,x,";1;");
+  hold on;
+  for i=2:5
+    [x,y] = hist(lsp(:,i) - lsp(:,i-1),100);
+    legend = sprintf(";%d;",i);
+    plot(y,x,legend);
+  endfor
+  hold off;
+  grid;
+
+  subplot(212)
+  [x,y] = hist(lsp(:,6)-lsp(:,5),100);
+  plot(y,x,";6;");
+  hold on;
+  for i=7:c
+    [x,y] = hist(lsp(:,i) - lsp(:,i-1),100);
+    legend = sprintf(";%d;",i);
+    plot(y,x,legend);
+  endfor
+  hold off;
+  grid;
+endfunction

diff --git a/codec2/octave/phase.m b/codec2/octave/phase.m
index fc13b47bf54a19c30e8ae2fd6b3ad142e05a1359..d9bab90646bb3a474b6c1bc3f80e73e0243f8651 100644 (file)
--- a/codec2/octave/phase.m
+++ b/codec2/octave/phase.m
@@ -2,24 +2,34 @@
 % David Rowe August 2009
 % experiments with phase for sinusoidal codecs
 
-function phase(samname, F0, randphase, png)
+function phase(samname, F0, png)
   Wo=2*pi*F0/8000;
   P=2*pi/Wo;
   L = floor(pi/Wo);
-  N = 16000;
+  Nsam = 16000;
+  N = 80;
+  F = Nsam/N;
   A = 10000/L;
   phi = zeros(1,L);
-  for m=1:L
-    phi(m) = m*Wo*50.125;
+  s = zeros(1,Nsam);
+  
+  for m=floor(L/2):L
+    phi_off(m) = -m*Wo*8;
   end
-  if (randphase == 1) 
-    rand("seed",0);
-    phi = rand(L,1)*2*pi;
-  endif
-  s = zeros(1,N);
 
-  for m=1:L
-    s = s + A*cos(m*Wo*(0:(N-1)) + phi(m));
+  for f=1:F
+    phi(1) = phi(1) + Wo*N;
+    phi(1) = mod(phi(1),2*pi);
+  
+    for m=2:L
+      phi(m) = m*phi(1);
+    end
+
+    x = zeros(1,N);
+    for m=1:L
+      x = x + A*cos(m*Wo*(0:(N-1)) + phi(m) + phi_off(m));
+    endfor
+    s((f-1)*N+1:f*N) = x;
   endfor
 
   figure(1);
@@ -30,7 +40,7 @@ function phase(samname, F0, randphase, png)
   fwrite(fs,s,"short");
   fclose(fs);
 
-  if png
+  if (nargin == 3)
       % small image to fit blog
 
       __gnuplot_set__ terminal png size 450,300

diff --git a/codec2/octave/phase2.m b/codec2/octave/phase2.m
new file mode 100644 (file)
index 0000000..5539bf1
--- /dev/null
+++ b/codec2/octave/phase2.m
@@ -0,0 +1,47 @@
+% phase2.m
+% David Rowe Sep 2009
+% experiments with phase for sinusoidal codecs, looking at phase
+% of excitation with real Am samples from hts1
+
+function phase2(samname, png)
+  N = 16000;
+
+  f=45;
+  model = load("../src/hts1a_model.txt");
+  Wo = model(f,1);
+  P=2*pi/Wo;
+  L = model(f,2);
+  A = model(f,3:(L+2));
+  
+  phi = zeros(1,L);
+  for m=1:L
+    phi(m) = -m*Wo*0.3*rand(1,1)*L;
+  end
+  s = zeros(1,N);
+
+  for m=1:L
+    s = s + A(m)*cos(m*Wo*(0:(N-1)) + phi(m));
+  endfor
+
+  figure(1);
+  clg;
+  plot(s(1:250));
+
+  fs=fopen(samname,"wb");
+  fwrite(fs,s,"short");
+  fclose(fs);
+
+  if (nargin == 2)
+      % small image to fit blog
+
+      __gnuplot_set__ terminal png size 450,300
+      ss = sprintf("__gnuplot_set__ output \"%s.png\"", samname);
+      eval(ss)
+      replot;
+
+      % for some reason I need this to stop large plot getting wiped
+      __gnuplot_set__ output "/dev/null"
+  endif
+
+endfunction
+

diff --git a/codec2/octave/plamp.m b/codec2/octave/plamp.m
index e4470699fc47018c490bce9a4f7346aea6729edb..44bd98af1425298fed1dd9de4fff6fcbde0d63d7 100644 (file)
--- a/codec2/octave/plamp.m
+++ b/codec2/octave/plamp.m
@@ -75,8 +75,8 @@ function plamp(samname, f)
       endif    
       signal = Am * Am';
       noise = (Am-Amq) * (Am-Amq)'; 
-      snr = 10*log10(signal/noise);
-      Am_err_label = sprintf(";Am_err SNR %4.2f dB;",snr);
+      snr1 = 10*log10(signal/noise);
+      Am_err_label = sprintf(";Am_err SNR %4.2f dB;",snr1);
       plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label);
     endif