From 688a2c278f60662c187a921a41fe2e25caffa07c Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Thu, 30 Aug 2012 07:37:15 +0000
Subject: [PATCH] various octave scripts develped during quality enhancement
 work over last few months

git-svn-id: https://svn.code.sf.net/p/freetel/code@650 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/octave/av_imp.m      |  43 ++++++++++
 codec2-dev/octave/cbphase.m     |  98 ++++++++++++++++++++++
 codec2-dev/octave/glottal.m     |   4 +
 codec2-dev/octave/lpcpf.m       |  46 +++++++++++
 codec2-dev/octave/phasesecord.m |  47 +++++++++++
 codec2-dev/octave/plamp.m       |  57 +++++++------
 codec2-dev/octave/pllpcpf.m     | 140 ++++++++++++++++++++++++++++++++
 codec2-dev/octave/plphase.m     |   2 +-
 codec2-dev/octave/plppe.m       |  65 +++++++++++++++
 codec2-dev/octave/sd.m          |  92 +++++++++++++++++++++
 codec2-dev/octave/twotone.m     |  52 ++++++++++++
 codec2-dev/octave/twotone1.m    |  76 +++++++++++++++++
 12 files changed, 695 insertions(+), 27 deletions(-)
 create mode 100644 codec2-dev/octave/av_imp.m
 create mode 100644 codec2-dev/octave/cbphase.m
 create mode 100644 codec2-dev/octave/lpcpf.m
 create mode 100644 codec2-dev/octave/phasesecord.m
 create mode 100644 codec2-dev/octave/pllpcpf.m
 create mode 100644 codec2-dev/octave/plppe.m
 create mode 100644 codec2-dev/octave/sd.m
 create mode 100644 codec2-dev/octave/twotone.m
 create mode 100644 codec2-dev/octave/twotone1.m

diff --git a/codec2-dev/octave/av_imp.m b/codec2-dev/octave/av_imp.m
new file mode 100644
index 00000000..8b7fa608
--- /dev/null
+++ b/codec2-dev/octave/av_imp.m
@@ -0,0 +1,43 @@
+% av_imp.m
+% David Rowe Aug 2012
+% Averages the impulse response samples
+
+function imp = av_imp(imp_filename, period_in_secs, st, en)
+  f = fopen(imp_filename,"rb");
+  s = fread(f, Inf, "short")';
+
+  Fs = 8000;
+  n = period_in_secs * Fs;
+
+  [r c] = size(s);
+
+  imp = zeros(1,n);
+  for i=1:n:c-n
+    imp = imp + s(i:i+n-1);
+  endfor
+  
+  % user supplies start and end samples after viweing plot
+
+  if (nargin == 4)
+    imp = imp(st:en);
+  end
+
+  % normalise
+
+  imp /= sqrt(sum(imp .^ 2));
+
+  [h w] = freqz(imp, 1, 4000);
+
+  figure(1);
+  clf;
+  plot(imp);
+
+  figure(2);
+  clf;
+  subplot(211)
+  plot(10*log10(abs(h)))
+  subplot(212)
+  plot(angle(h))
+
+endfunction
+
diff --git a/codec2-dev/octave/cbphase.m b/codec2-dev/octave/cbphase.m
new file mode 100644
index 00000000..8e82da1c
--- /dev/null
+++ b/codec2-dev/octave/cbphase.m
@@ -0,0 +1,98 @@
+% cbphase.m
+% David Rowe Aug 2012
+% Used to experiment with critical band phase perception and smoothing
+
+function cbphase
+
+  Wo = 100.0*pi/4000;
+  L = floor(pi/Wo);
+
+  A = zeros(1,L);
+  phi = zeros(1,L);
+
+  % three harmonics in this band
+
+  b = 4; a = b-1; c = b+1;
+
+  % set up phases and mags for 2nd order system (see phasesecord.m)
+   
+  wres = b*Wo;
+  phi(a) = 3*pi/4 + wres;
+  phi(b) = pi/2 + wres;
+  phi(c) = pi/4 + wres;
+
+  A(a) = 0.707;
+  A(b) = 1;
+  A(c) = 0.707;
+
+  % add linear component
+
+  phi(1) = pi;
+  phi(2:L) = phi(2:L) + (2:L)*phi(1);
+  phi = phi - 2*pi*(floor(phi/(2*pi)) + 0.5);
+
+  N = 16000;
+  Nplot = 250;
+  s = zeros(1,N);
+
+  for m=a:c
+    s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi(m));
+    s = s + s_m;
+  endfor
+
+  figure(2);
+  clf;
+  subplot(211)
+  plot((1:L)*Wo*4000/pi, A,'+');
+  subplot(212)
+  plot((1:L)*Wo*4000/pi, phi,'+');
+
+  %v = A(a)*exp(j*phi(a)) + A(b)*exp(j*phi(b)) + A(c)*exp(j*phi(c));
+  %compass(v,"r")
+  %hold off;
+  
+  % est phi1
+
+  diff = phi(b) - phi(a)
+  sumi = sin(diff);
+  sumr = cos(diff);
+  diff = phi(c) - phi(b)
+  sumi += sin(diff);
+  sumr += cos(diff);
+  phi1_ = atan2(sumi, sumr)
+  s_v = cos(Wo*(0:(N-1)) + phi1_);
+
+  figure(1);
+  clf;
+  subplot(211)
+  plot(s(1:Nplot));
+  hold on;
+  plot(s_v(1:Nplot),"r");
+  hold off;
+
+  % build (hopefully) perceptually similar phase
+
+  phi_(a) = a*phi1_;
+  phi_(b) = b*phi1_;
+  phi_(c) = c*phi1_;
+
+  s_ = zeros(1,N);
+
+  for m=a:c
+    s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi_(m));
+    s_ = s_ + s_m;
+  endfor
+ 
+  subplot(212)
+  plot(s_(1:Nplot));
+  
+  gain = 8000;
+  fs=fopen("orig_ph.raw","wb");
+  fwrite(fs,gain*s,"short");
+  fclose(fs);
+  fs=fopen("mod_ph.raw","wb");
+  fwrite(fs,gain*s_,"short");
+  fclose(fs);
+
+endfunction
+
diff --git a/codec2-dev/octave/glottal.m b/codec2-dev/octave/glottal.m
index 23f8df97..46675e7d 100644
--- a/codec2-dev/octave/glottal.m
+++ b/codec2-dev/octave/glottal.m
@@ -11,10 +11,14 @@ kexc = [ 8,  -16,   26, -48,  86, -162, 294, -502, 718, -728, 184 672, -610, -67
 kexc = shift(kexc,sh);
 kexc = [kexc(1:sh) zeros(1,512-25) kexc(sh+1:25)];
 figure(1)
+clf
 plot(kexc)
 figure(2)
 G = fft(kexc);
+subplot(211)
 plot((1:256)*(4000/256),unwrap(angle(G(1:256))))
+subplot(212)
+plot(20*log10(abs(G)))
 
 f=fopen("glottal.c","wt");
 fprintf(f,"const float glottal[]={\n");
diff --git a/codec2-dev/octave/lpcpf.m b/codec2-dev/octave/lpcpf.m
new file mode 100644
index 00000000..f1e0982c
--- /dev/null
+++ b/codec2-dev/octave/lpcpf.m
@@ -0,0 +1,46 @@
+% lpcpf.m
+% David Rowe Aug 27 2012
+% Experiments with LPC post filtering
+
+function lpcpf(ak_filename, f)
+  aks = load(ak_filename);
+	 
+  ak = aks(f,:);
+  [tmp p] = size(ak);
+  p -= 1;
+
+  A = freqz(1,ak, 4000);	 
+  AdB = 20*log10(abs(A));
+
+  gamma = 0.5;
+  gammas = gamma .^ (0:p);
+  W = freqz(ak .* gammas,1, 4000);
+  WdB = 20*log10(abs(W));
+
+  beta = 0.2;
+  R = abs(freqz(ak .* gammas, ak, 4000));
+  %P = (R/max(R)) .^ beta;
+  P = R .^ beta;
+  AP = abs(A) .* P;
+
+  eA = sum(abs(A) .^ 2);
+  eAP = sum(AP .^ 2);
+  gain = sqrt(eA/eAP)
+  AP *= gain;
+
+  PdB = 20*log10(P);
+
+  APdB = 20*log10(AP);
+  10*log10(sum(AP .^ 2))/10*log10(sum(abs(A) .^ 2))
+
+  figure(1);
+  clf;
+  plot(AdB);
+  hold on;
+  plot(WdB,'g');
+  plot(PdB,'r');
+  plot(APdB,'b.');
+  hold off;
+
+endfunction
+
diff --git a/codec2-dev/octave/phasesecord.m b/codec2-dev/octave/phasesecord.m
new file mode 100644
index 00000000..a3cf2516
--- /dev/null
+++ b/codec2-dev/octave/phasesecord.m
@@ -0,0 +1,47 @@
+% phasesecord.m
+% David Rowe Aug 2012
+% Used to experiment with aproximations of phase of 2nd order systems
+
+function phasesecord(w,beta)
+
+  a = [1 -2*cos(w)*beta beta*beta];
+  b = 1;
+
+  [h w1] = freqz(b,a);
+
+  figure(1)
+  subplot(211)
+  plot(abs(h))
+  subplot(212)
+  plot(angle(h))
+
+  % for beta close to 1, we approximate 3 dB points as 1-beta above
+  % and below the resonance freq. Note this fails if w=0 as there is a
+  % double pole.  Lets sample the freq response at the 3dB points and
+  % w:
+
+  ws = [w-(1-beta) w w+(1-beta)];
+  [h w1] = freqz(b,a,ws);
+
+  % gain as a fraction of max, should be 3dB. Within 1.3 dB or for w > pi/8,
+  % gets innacurate near w=0 due to 2nd pole
+
+  printf("mag measured...:"); printf("% 4.3f ", abs(h)/max(abs(h)));
+
+  % measured angle, 45 deg from angle at w
+
+  printf("\nangle measured.: "); printf("% 5.3f ", angle(h));
+
+  % Our estimate of angle, (pi+w) is phase at resonance, at lower 3dB
+  % phase is pi/4 ahead, at upper 3B pi/4 behind.  -pi/2 is contribution of
+  % other pole at at -w to phase
+
+  ph_lower = (pi+w) + pi/4 - pi/2;
+  ph_res   =(pi+w) - pi/2;
+  ph_upper = (pi+w) - pi/4 - pi/2;
+  ph_ests = [ph_lower ph_res ph_upper];
+  ph_ests = ph_ests - 2*pi*(floor(ph_ests/(2*pi)) + 0.5);
+  printf("\nangle estimated:"); printf("% 5.3f ", ph_ests);
+  printf("\n");
+endfunction
+
diff --git a/codec2-dev/octave/plamp.m b/codec2-dev/octave/plamp.m
index 79683596..62b6893a 100644
--- a/codec2-dev/octave/plamp.m
+++ b/codec2-dev/octave/plamp.m
@@ -4,8 +4,16 @@
 %
 % Plot ampltiude modelling information from dump files.
 
-function plamp(samname, f)
+function plamp(samname, f, samname2)
   
+  % switch some stuff off to unclutter display
+
+  plot_lsp = 0;
+  plot_snr = 0;
+  plot_vsnr = 0;
+  plot_sw = 0;
+  plot_pw = 0;
+
   sn_name = strcat(samname,"_sn.txt");
   Sn = load(sn_name);
 
@@ -60,6 +68,16 @@ function plamp(samname, f)
     snr = load(snr_name);
   endif
 
+  % optional second file, for exploring post filter
+
+  model2q_name = " ";
+  if nargin == 3
+    model2q_name = strcat(samname2,"_qmodel.txt");
+    if file_in_path(".",modelq_name)
+      model2q = load(model2q_name);
+    end
+  end
+
   Ew_on = 1;
   k = ' ';
   do 
@@ -67,7 +85,7 @@ function plamp(samname, f)
     clf;
 %    s = [ Sn(2*(f-2)-1,:) Sn(2*(f-2),:) ];
     s = [ Sn(2*f-1,:) Sn(2*f,:) ];
-    size(s)
+    size(s);
     plot(s);
     axis([1 length(s) -20000 20000]);
 
@@ -78,20 +96,14 @@ function plamp(samname, f)
     plot((1:L)*Wo*4000/pi, 20*log10(Am),";Am;r");
     axis([1 4000 -10 80]);
     hold on;
-%    plot((0:255)*4000/256, Sw(f-2,:),";Sw;");
-    plot((0:255)*4000/256, Sw(f,:),";Sw;");
-    
-    if (file_in_path(".",sw__name))
-       plot((0:255)*4000/256, Sw_(f,:),";Sw_;g");
-       if (Ew_on == 1)
-          plot((0:255)*4000/256, Ew(f,:),";Ew_;r");
-       endif
-    endif
+    if plot_sw
+      plot((0:255)*4000/256, Sw(f,:),";Sw;");
+    end
  
     if (file_in_path(".",modelq_name))
       Amq = modelq(f,3:(L+2));
       plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;g" );
-      if (file_in_path(".",pw_name))
+      if (file_in_path(".",pw_name) && plot_pw)
         plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;c");
       endif	
       signal = Am * Am';
@@ -101,7 +113,12 @@ function plamp(samname, f)
       plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label);
     endif
 
-    if (file_in_path(".",snr_name))
+    if file_in_path(".",model2q_name)
+      Amq2 = model2q(f,3:(L+2));
+      plot((1:L)*Wo*4000/pi, 20*log10(Amq2),";Amq2;m" );
+    end
+
+    if (file_in_path(".",snr_name) && plot_vsnr)
       snr_label = sprintf(";Voicing SNR %4.2f dB;",snr(f));
       plot(1,1,snr_label);
     endif
@@ -119,7 +136,7 @@ function plamp(samname, f)
       %plot((1:L)*Wo*4000/pi, 20*log10(orig-synth), phase_err_label);
     endif
 
-    if (file_in_path(".",lsp_name))
+    if (file_in_path(".",lsp_name) && plot_lsp)
       for l=1:10
         plot([lsp(f,l)*4000/pi lsp(f,l)*4000/pi], [60 80], 'r');
       endfor
@@ -139,18 +156,6 @@ function plamp(samname, f)
       %figure(2);
     %endif
 
-    % autocorrelation function to research voicing est
-    
-    if (file_in_path(".",rk_name))
-      figure(3);
-      plot(Rk(f,:) / Rk(f,1), ";Rk;");
-      hold on;
-      p = floor(2*pi/Wo);
-      plot([p p ], [0 1], 'r');
-      hold off;
-      %figure(2);
-    endif
-
     % interactive menu
 
     printf("\rframe: %d  menu: n-next  b-back  p-png  q-quit e-toggle Ew", f);
diff --git a/codec2-dev/octave/pllpcpf.m b/codec2-dev/octave/pllpcpf.m
new file mode 100644
index 00000000..b3912f8a
--- /dev/null
+++ b/codec2-dev/octave/pllpcpf.m
@@ -0,0 +1,140 @@
+% Copyright David Rowe 2012
+% This program is distributed under the terms of the GNU General Public License 
+% Version 2
+%
+% Plot amplitude modelling information from dump files to test and develop
+% LPC post filter.
+
+function pllpcpf(samname, f)
+  
+  % switch some stuff off to unclutter display
+
+  plot_lsp = 0;
+  plot_snr = 0;
+  plot_vsnr = 0;
+  plot_sw = 1;
+  plot_pw = 1;
+  plot_pwb = 1;
+  plot_rw = 1;
+
+  sn_name = strcat(samname,"_sn.txt");
+  Sn = load(sn_name);
+
+  sw_name = strcat(samname,"_sw.txt");
+  Sw = load(sw_name);
+
+  sw__name = strcat(samname,"_sw_.txt");
+  if (file_in_path(".",sw__name))
+    Sw_ = load(sw__name);
+  endif
+
+  model_name = strcat(samname,"_model.txt");
+  model = load(model_name);
+
+  modelq_name = strcat(samname,"_qmodel.txt");
+  if (file_in_path(".",modelq_name))
+    modelq = load(modelq_name);
+  endif
+
+  % Pw (LPC synth filter spectrum) before post filter
+
+  pwb_name = strcat(samname,"_pwb.txt");
+  if (file_in_path(".",pwb_name))
+    Pwb = load(pwb_name);
+  endif
+
+  % Rw (Post filter spectrum)
+
+  rw_name = strcat(samname,"_rw.txt");
+  if (file_in_path(".",rw_name))
+    Rw = load(rw_name);
+  endif
+
+  % Pw (LPC synth filter spectrum) after post filter
+
+  pw_name = strcat(samname,"_pw.txt");
+  if (file_in_path(".",pw_name))
+    Pw = load(pw_name);
+  endif
+
+
+  Ew_on = 1;
+  k = ' ';
+  do 
+    figure(1);
+    clf;
+    s = [ Sn(2*f-1,:) Sn(2*f,:) ];
+    size(s);
+    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;r");
+    axis([1 4000 -10 80]);
+    hold on;
+    if plot_sw
+      plot((0:255)*4000/256, Sw(f,:),";Sw;");
+    end
+ 
+    if (file_in_path(".",modelq_name))
+
+      Amq = modelq(f,3:(L+2));
+      plot((1:L)*Wo*4000/pi, 20*log10(Amq),";Amq;g" );
+
+      if (file_in_path(".",pwb_name) && plot_pwb)
+        plot((0:255)*4000/256, 10*log10(Pwb(f,:)),";Pwb;c");
+      endif	
+
+      if (file_in_path(".",rw_name) && plot_rw)
+        plot((0:255)*4000/256, 10*log10(Rw(f,:)),";Rw;c");
+      endif	
+
+      if (file_in_path(".",pw_name) && plot_pw)
+        plot((0:255)*4000/256, 10*log10(Pw(f,:)),";Pw;c.");
+      endif	
+
+      signal = Am * Am';
+      noise = (Am-Amq) * (Am-Amq)'; 
+      snr1 = 10*log10(signal/noise);
+      Am_err_label = sprintf(";Am error SNR %4.2f dB;m",snr1);
+      plot((1:L)*Wo*4000/pi, 20*log10(Amq) - 20*log10(Am), Am_err_label);
+    endif
+
+
+    hold off;
+
+    % interactive menu
+
+    printf("\rframe: %d  menu: n-next  b-back  p-png  q-quit", f);
+    fflush(stdout);
+    k = kbhit();
+    if (k == 'n')
+      f = f + 1;
+    endif
+    if (k == 'b')
+      f = f - 1;
+    endif
+
+    % optional print to PNG
+
+    if (k == 'p')
+      figure(1);
+      pngname = sprintf("%s_%d_sn.png",samname,f);
+      print(pngname, '-dpng', "-S500,500")
+      pngname = sprintf("%s_%d_sn_large.png",samname,f);
+      print(pngname, '-dpng', "-S800,600")
+
+      figure(2);
+      pngname = sprintf("%s_%d_sw.png",samname,f);
+      print(pngname, '-dpng', "-S500,500")
+      pngname = sprintf("%s_%d_sw_large.png",samname,f);
+      print(pngname, '-dpng', "-S1200,800")
+    endif
+
+  until (k == 'q')
+  printf("\n");
+
+endfunction
diff --git a/codec2-dev/octave/plphase.m b/codec2-dev/octave/plphase.m
index 9e611856..c12422ea 100644
--- a/codec2-dev/octave/plphase.m
+++ b/codec2-dev/octave/plphase.m
@@ -112,7 +112,7 @@ function plphase(samname, f)
       axis;
       if (file_in_path(".", phase_name_))
         hold on;
-        plot((1:L)*Wo*4000/pi, phase_(f,1:L)*180/pi, "g;phase_;");
+        plot((1:L)*Wo*4000/pi, phase_(f,1:L)*180/pi, "g;phase after;");
 	grid
 	hold off;
       endif
diff --git a/codec2-dev/octave/plppe.m b/codec2-dev/octave/plppe.m
new file mode 100644
index 00000000..cbc5b562
--- /dev/null
+++ b/codec2-dev/octave/plppe.m
@@ -0,0 +1,65 @@
+% Copyright David Rowe 2009
+% This program is distributed under the terms of the GNU General Public License 
+% Version 2
+%
+% Plot two sparse phase prediction error text files.
+% Generate data from print_pred_error, print_pred_error_sparse_wo_correction1 etc
+
+function plppe(ppe1_file, ppe2_file, f)
+  
+  ppe1 = load(ppe1_file);
+  ppe2 = load(ppe2_file);
+
+  std1 = std(nonzeros(ppe1(:,40:80)));
+  std2 = std(nonzeros(ppe2(:,40:80)));
+
+  printf("std dev for %s is %4.3f\n", ppe1_file, std1);
+  printf("std dev for %s is %4.3f\n", ppe2_file, std2);
+ 
+  figure(1);
+  clf;
+  subplot(211)
+  hist(nonzeros(ppe1(:,40:80)),20);
+  subplot(212)
+  hist(nonzeros(ppe2(:,40:80)),20);
+  
+  k = ' ';
+  do 
+    figure(2);
+    clf;
+    subplot(211)
+    L = length(nonzeros(ppe1(f,:)));
+    x = (1:L)*4000/L;
+    std1 = std(nonzeros(ppe1(f,:)));
+    legend = sprintf(";std dev %4.3f;", std1);
+    plot(x, nonzeros(ppe1(f,:)),legend);
+    axis([0 4000 -pi pi]);
+    subplot(212)
+    std2 = std(nonzeros(ppe2(f,:)));
+    legend = sprintf(";std dev %4.3f;", std2);
+    plot(x, nonzeros(ppe2(f,:)),legend);
+    axis([0 4000 -pi pi]);
+
+    % interactive menu
+
+    printf("\rframe: %d  menu: n-next  b-back  p-png  q-quit ", f);
+    fflush(stdout);
+    k = kbhit();
+    if (k == 'n')
+      f = f + 1;
+    endif
+    if (k == 'b')
+      f = f - 1;
+    endif
+
+    % optional print to PNG
+
+    if (k == 'p')
+       pngname = sprintf("%s_%d",samname,f);
+       png(pngname);
+    endif
+
+  until (k == 'q')
+  printf("\n");
+
+endfunction
diff --git a/codec2-dev/octave/sd.m b/codec2-dev/octave/sd.m
new file mode 100644
index 00000000..112af614
--- /dev/null
+++ b/codec2-dev/octave/sd.m
@@ -0,0 +1,92 @@
+% sd.m
+% David Rowe Aug 2012
+% Plots the spectal distorion between twofiles of LPCs.  Used for LSP 
+% quantisation tuning.
+
+function sd(raw_filename, dump_file_prefix, f)
+
+  ak1_filename = sprintf("%s_ak.txt", dump_file_prefix);
+  ak2_filename = sprintf("%s_ak_.txt", dump_file_prefix);
+  ak1 = load(ak1_filename);
+  ak2 = load(ak2_filename);
+
+  [ak1_r, ak1_c] = size(ak1);
+  [ak2_r, ak2_c] = size(ak1);
+
+  frames = max([ak1_r ak2_r]);
+  sd = zeros(1,frames);
+  Ndft = 512;
+  A1 = zeros(frames, Ndft);
+  A2 = zeros(frames, Ndft);
+
+  % initial helicopter view of all frames
+
+  for i = 1:frames
+    A1(i,:) = -20*log10(abs(fft(ak1(i,:),Ndft)));
+    A2(i,:) = -20*log10(abs(fft(ak2(i,:),Ndft)));
+    sd(i) = sum((A1(i,:) - A2(i,:)).^2)/Ndft;
+  end
+
+  figure(1);
+  clf;
+  subplot(211)
+  fs=fopen(raw_filename,"rb");
+  s = fread(fs,Inf,"short");
+  plot(s);
+  subplot(212)
+  plot(sd);
+
+  % now enter single step mode so we can analyse each frame
+
+  sn_name = strcat(dump_file_prefix,"_sn.txt");
+  Sn = load(sn_name);
+
+  lsp1_filename = sprintf("%s_lsp.txt", dump_file_prefix);
+  lsp2_filename = sprintf("%s_lsp_.txt", dump_file_prefix);
+  lsp1 = load(lsp1_filename);
+  lsp2 = load(lsp2_filename);
+
+  k = ' ';
+  do 
+    figure(2);
+    clf;
+    subplot(211)
+    s = [ Sn(2*f-1,:) Sn(2*f,:) ];
+    size(s);
+    plot(s);
+    axis([1 length(s) -20000 20000]);
+
+    subplot(212);
+    plot((1:Ndft/2)*4000/(Ndft/2), A1(f,1:(Ndft/2)),";A1;r");
+    axis([1 4000 -20 40]);
+    hold on;
+    plot((1:Ndft/2)*4000/(Ndft/2), A2(f,1:(Ndft/2)),";A2;");
+ 
+    for l=1:10
+        plot([lsp1(f,l)*4000/pi lsp1(f,l)*4000/pi], [0  -10], 'r');
+        plot([lsp2(f,l)*4000/pi lsp2(f,l)*4000/pi], [-10 -20], 'b');
+    endfor
+    plot(0,0,';lsp1;r');
+    plot(0,0,';lsp2;b');
+    sd_str = sprintf(";sd %3.2f dB*dB;", sd(f));
+    plot(0,0,sd_str);
+   
+    hold off;
+
+    % interactive menu
+
+    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");
+
+endfunction
+
diff --git a/codec2-dev/octave/twotone.m b/codec2-dev/octave/twotone.m
new file mode 100644
index 00000000..a20c5c5a
--- /dev/null
+++ b/codec2-dev/octave/twotone.m
@@ -0,0 +1,52 @@
+% twotone.m
+% David Rowe Aug 2012
+% Used to experiment with combining phase of two tones
+
+function cbphase
+
+  Wo = 100.0*pi/4000;
+  L = floor(pi/Wo);
+  phi = zeros(1,L);
+
+  % two harmonics 
+
+  a = 20; b = 21;
+
+  % set up phases to whatever
+   
+  phi(a) = -pi;
+  phi(b) = -pi/2;
+
+  % synthesis the two-tone signal
+
+  N = 16000;
+  Nplot = 250;
+  s = zeros(1,N);
+
+  for m=a:b
+    s_m = cos(m*Wo*(0:(N-1)) + phi(m));
+    s = s + s_m;
+  endfor
+
+  % now our theory says that this signal should be the same perceptually
+
+  phi_(a) = (phi(a) - phi(b))/2;
+  phi_(b) = (phi(b) - phi(a))/2;
+
+  s_ = zeros(1,N);
+  for m=a:b
+    s_m = cos(m*Wo*(0:(N-1)) + phi_(m));
+    s_ = s_ + s_m;
+  endfor
+
+  % plot them and see if envelope has the same phase, but "carriers"
+  % have different phase
+
+  figure(1);
+  clf;
+  subplot(211);
+  plot(s(1:Nplot));
+  subplot(212);
+  plot(s_(1:Nplot),'r');
+endfunction
+
diff --git a/codec2-dev/octave/twotone1.m b/codec2-dev/octave/twotone1.m
new file mode 100644
index 00000000..9f36c68d
--- /dev/null
+++ b/codec2-dev/octave/twotone1.m
@@ -0,0 +1,76 @@
+% twotone1.m
+% David Rowe Aug 17 2012
+%
+% Used to experiment with combining phase of two tones.  This version
+% sets up a complete synthetic speech signal then tries to combine the
+% phase of high frequency tones.  Lets see if we can do this and keep perceptual 
+% nature of signal the same.
+
+function twotone1
+ 
+  % hts1a frame 47
+
+  Wo = 0.093168;
+  L = 33;
+  %A = [69.626907 460.218536 839.677429 2577.498047 972.647888 712.755066 489.048553 364.830536 409.230652 371.767487 489.112854	893.127014 2447.596680 752.878113 475.720520 234.452271	248.161606 232.171051 202.669891 323.914490 678.749451 362.958038 211.652512 170.764435	148.631790 169.261673 272.254150 176.872375 67.344391 99.022301	60.812035 34.319073 14.864757];
+  A  = zeros(1,L)*100;
+  phi = [1.560274 1.508063 -1.565184 1.289117 -2.547365	1.412528 -1.303992 3.121130 1.087573 -1.158161 -2.928007 0.995093 -2.614023 0.246136 -2.267406 2.143802	-0.273431 -2.266897 1.685171 -0.668712 2.699722	-1.151891 2.406379 -0.046192 -2.718611 0.761067	-2.305014 0.133172 -1.428978 1.492630 -1.668385	1.539734 -1.336615];
+  %phi = zeros(1,L);
+  st = floor(L/2);
+  %st = 1;
+
+  A(st:st+5) = 1000;
+
+  % now set up phase of signal with phase of upper frequency harmonic
+  % pairs combined
+
+  phi_ = phi;
+  for m=floor(L/2):2:L
+    phi_(m)   = (phi(m) - phi(m+1))/2;
+    phi_(m+1) = (phi(m+1) - phi(m))/2;
+    %phi_(m+1) = 0;
+  end
+
+  % synthesise the signals
+
+  N = 16000;
+  Nplot = 250;
+
+  s = zeros(1,N);
+  for m=st:L
+    s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi(m));
+    s = s + s_m;
+  endfor
+
+  s_ = zeros(1,N);
+  for m=st:L
+    s_m = A(m)*cos(m*Wo*(0:(N-1)) + phi_(m));
+    s_ = s_ + s_m;
+  endfor
+
+  % plot them, expect to see similar time domain waveforms
+
+  figure(1);
+  clf;
+  subplot(211);
+  plot(s(1:Nplot));
+  subplot(212);
+  plot(s_(1:Nplot),'r');
+
+  figure(2);
+  clf;
+  subplot(211);
+  plot(s(1:Nplot)-s_(1:Nplot));
+
+  % save to disk
+
+  gain = 1;
+  fs=fopen("twotone1_orig.raw","wb");
+  fwrite(fs,gain*s,"short");
+  fclose(fs);
+  fs=fopen("twotone1_comb.raw","wb");
+  fwrite(fs,gain*s_,"short");
+  fclose(fs);
+  
+endfunction
+
-- 
2.25.1