From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Wed, 26 Aug 2009 05:00:03 +0000 (+0000)
Subject: lpc modelling coded, plamp.m script partially coded, menu.sh
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=764a2a4a72c66b0b2c653e2be46d55cc454d01ac;p=freetel-svn-tracking.git

lpc modelling coded, plamp.m script partially coded, menu.sh

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

diff --git a/codec2/octave/plamp.m b/codec2/octave/plamp.m
new file mode 100644
index 00000000..8ae51c6c
--- /dev/null
+++ b/codec2/octave/plamp.m
@@ -0,0 +1,34 @@
+% Copyright David Rowe 2009
+% This program is distributed under the terms of the GNU General Public License 
+% Version 2
+
+function plamp(samname, f)
+  
+  sn_name = strcat(samname,"_sn.txt");
+  Sn = load(sn_name);
+  model_name = strcat(samname,"_model.txt");
+  model = load(model_name);
+  modelq_name = strcat(samname,"_qmodel.txt");
+  modelq = load(modelq_name);
+  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;
+
+  Am_err = mean(abs(20*log10(Amq) - 20*log10(Am)))
+endfunction
diff --git a/codec2/pitch/hts2.p b/codec2/pitch/hts2.p
new file mode 100755
index 00000000..20e26804
--- /dev/null
+++ b/codec2/pitch/hts2.p
@@ -0,0 +1,300 @@
+  0.0000000e+000
+  9.2753623e+001
+  5.4237288e+001
+  8.5906040e+001
+  7.0329670e+001
+  5.5652174e+001
+  5.4237288e+001
+  5.4935622e+001
+  5.4700855e+001
+  7.5739645e+001
+  7.3563218e+001
+  1.2307692e+002
+  1.1428571e+002
+  7.3563218e+001
+  7.7108434e+001
+  1.8550725e+002
+  1.2673267e+002
+  1.0847458e+002
+  7.8527607e+001
+  8.8888889e+001
+  8.3116883e+001
+  8.1012658e+001
+  1.0756303e+002
+  1.3061224e+002
+  4.8301887e+001
+  4.7940075e+001
+  4.8120301e+001
+  4.9230769e+001
+  4.9420849e+001
+  4.6886447e+001
+  4.2953020e+001
+  3.9263804e+001
+  3.7869822e+001
+  3.5457064e+001
+  3.4224599e+001
+  3.3333333e+001
+  3.2820513e+001
+  3.2000000e+001
+  3.1295844e+001
+  2.9906542e+001
+  2.9493088e+001
+  2.9090909e+001
+  2.8699552e+001
+  2.8131868e+001
+  2.7826087e+001
+  2.7826087e+001
+  2.7826087e+001
+  2.8193833e+001
+  2.7467811e+001
+  2.6890756e+001
+  5.4468085e+001
+  5.4237288e+001
+  6.4974619e+001
+  1.0756303e+002
+  8.8888889e+001
+  1.0406504e+002
+  4.4599303e+001
+  5.4468085e+001
+  3.6260623e+001
+  3.6260623e+001
+  8.1012658e+001
+  7.0329670e+001
+  1.2929293e+002
+  9.9224806e+001
+  4.3097643e+001
+  4.4137931e+001
+  4.5714286e+001
+  4.7407407e+001
+  4.8301887e+001
+  4.9230769e+001
+  4.9420849e+001
+  5.0996016e+001
+  5.1405622e+001
+  5.1405622e+001
+  5.2244898e+001
+  5.2459016e+001
+  5.2459016e+001
+  5.2244898e+001
+  5.3333333e+001
+  5.2459016e+001
+  5.2244898e+001
+  5.1405622e+001
+  5.1405622e+001
+  5.1200000e+001
+  5.0996016e+001
+  5.0196078e+001
+  4.9230769e+001
+  4.9230769e+001
+  4.9230769e+001
+  4.9420849e+001
+  4.9230769e+001
+  4.9042146e+001
+  9.8461538e+001
+  1.0158730e+002
+  5.1821862e+001
+  9.0140845e+001
+  1.0491803e+002
+  1.4382022e+002
+  5.2459016e+001
+  5.2459016e+001
+  1.2929293e+002
+  1.6410256e+002
+  8.0000000e+001
+  7.3563218e+001
+  1.0158730e+002
+  9.9224806e+001
+  4.9042146e+001
+  4.9042146e+001
+  4.9042146e+001
+  5.9259259e+001
+  1.4382022e+002
+  7.2316384e+001
+  1.0847458e+002
+  1.1228070e+002
+  1.6202532e+002
+  8.1528662e+001
+  7.2727273e+001
+  1.8550725e+002
+  6.0093897e+001
+  1.0847458e+002
+  8.9510490e+001
+  7.1508380e+001
+  4.0125392e+001
+  4.0634921e+001
+  4.0634921e+001
+  4.0251572e+001
+  4.0506329e+001
+  4.3986254e+001
+  4.0506329e+001
+  9.8461538e+001
+  5.6140351e+001
+  6.5641026e+001
+  5.4237288e+001
+  1.1636364e+002
+  3.4316354e+001
+  3.4972678e+001
+  3.7758112e+001
+  4.0634921e+001
+  4.0506329e+001
+  4.1290323e+001
+  4.2524917e+001
+  4.3389831e+001
+  4.4599303e+001
+  4.4912281e+001
+  4.6545455e+001
+  4.7232472e+001
+  4.8301887e+001
+  4.9230769e+001
+  4.9420849e+001
+  5.0393701e+001
+  5.1405622e+001
+  5.3333333e+001
+  5.3112033e+001
+  1.1034483e+002
+  9.7709924e+001
+  1.4382022e+002
+  5.0996016e+001
+  5.1821862e+001
+  5.0996016e+001
+  5.2032520e+001
+  5.3112033e+001
+  5.3556485e+001
+  5.4468085e+001
+  5.5652174e+001
+  5.4700855e+001
+  5.4700855e+001
+  5.4935622e+001
+  5.4700855e+001
+  5.4700855e+001
+  5.4468085e+001
+  5.4468085e+001
+  5.4468085e+001
+  5.4468085e+001
+  5.3333333e+001
+  5.1405622e+001
+  5.0996016e+001
+  5.0000000e+001
+  4.8120301e+001
+  4.8669202e+001
+  4.7058824e+001
+  4.6376812e+001
+  4.5070423e+001
+  4.4912281e+001
+  4.4137931e+001
+  4.2809365e+001
+  4.2666667e+001
+  4.2105263e+001
+  4.1423948e+001
+  4.1290323e+001
+  4.1290323e+001
+  4.1290323e+001
+  4.0634921e+001
+  4.0634921e+001
+  4.0634921e+001
+  4.0634921e+001
+  4.0764331e+001
+  4.1423948e+001
+  4.2953020e+001
+  4.5551601e+001
+  1.7534247e+002
+  4.7232472e+001
+  1.3763441e+002
+  1.3061224e+002
+  4.5551601e+001
+  4.3686007e+001
+  4.8669202e+001
+  9.4117647e+001
+  8.1012658e+001
+  1.1228070e+002
+  1.3617021e+002
+  4.3097643e+001
+  4.3835616e+001
+  4.6376812e+001
+  4.6545455e+001
+  4.6043165e+001
+  4.8301887e+001
+  4.9042146e+001
+  4.9420849e+001
+  5.1200000e+001
+  5.1405622e+001
+  5.2244898e+001
+  1.2929293e+002
+  1.2929293e+002
+  1.5238095e+002
+  1.5238095e+002
+  1.3913043e+002
+  9.0140845e+001
+  1.0940171e+002
+  9.0140845e+001
+  1.2307692e+002
+  8.9510490e+001
+  6.9565217e+001
+  7.3142857e+001
+  1.1034483e+002
+  7.8048780e+001
+  7.2727273e+001
+  1.0078740e+002
+  1.0940171e+002
+  1.1743119e+002
+  8.7074830e+001
+  1.8550725e+002
+  6.5306122e+001
+  1.3617021e+002
+  5.2674897e+001
+  1.0940171e+002
+  1.5238095e+002
+  1.4065934e+002
+  1.0756303e+002
+  1.0406504e+002
+  5.0793651e+001
+  4.9420849e+001
+  4.4444444e+001
+  7.0329670e+001
+  7.2727273e+001
+  7.4418605e+001
+  1.1636364e+002
+  1.0406504e+002
+  1.2307692e+002
+  1.2549020e+002
+  1.7297297e+002
+  4.5878136e+001
+  4.9805447e+001
+  6.2745098e+001
+  9.2086331e+001
+  9.1428571e+001
+  5.7142857e+001
+  4.8484848e+001
+  4.1157556e+001
+  2.2857143e+001
+  3.0046948e+001
+  9.4814815e+001
+  5.7918552e+001
+  9.0140845e+001
+  7.4418605e+001
+  7.4418605e+001
+  5.4700855e+001
+  9.5522388e+001
+  7.4853801e+001
+  9.4117647e+001
+  9.5522388e+001
+  9.9224806e+001
+  8.1012658e+001
+  1.1851852e+002
+  6.8817204e+001
+  8.5906040e+001
+  6.7015707e+001
+  4.3537415e+001
+  6.5306122e+001
+  3.1295844e+001
+  7.5739645e+001
+  6.2135922e+001
+  9.9224806e+001
+  5.7657658e+001
+  5.2244898e+001
+  5.8447489e+001
+  0.0000000e+000
+  0.0000000e+000
+  0.0000000e+000
+  0.0000000e+000
+  0.0000000e+000
diff --git a/codec2/script/menu.sh b/codec2/script/menu.sh
new file mode 100755
index 00000000..ce9a89d9
--- /dev/null
+++ b/codec2/script/menu.sh
@@ -0,0 +1,66 @@
+#!/bin/sh
+# ./menu.sh
+#
+# David Rowe
+# Created August 2009
+#
+# Presents a menu of sound files, press 1 to play file1, 2 to play file2 etc
+#
+# The aim is to make comparing files with different processing easier than
+# using up-arrow on the command line.  Based on cdialog.
+#
+# usage:
+#   menu.sh file1.raw file2.raw ........ [-d playbackdevice]
+#
+# for example:
+#
+#   ../script/menu.sh hts1a.raw hts1a_uq.raw 
+#
+# or:
+#
+#   ../script/menu.sh hts1a.raw hts1a_uq.raw -d /dev/dsp1
+#
+
+#  Copyright (C) 2007 David Rowe
+# 
+#  All rights reserved.
+# 
+#  This program is free software; you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License version 2, as
+#  published by the Free Software Foundation.
+# 
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+# 
+#  You should have received a copy of the GNU General Public License
+#  along with this program; if not, write to the Free Software
+#  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+files=0
+items="0 QUIT "
+while [ ! -z "$1" ]
+do
+  case "$1" in
+    -d) dsp="${1} ${2}"; shift;;
+     *) files=`expr 1 + $files`;
+        new_file=$1;
+        file[$files]=$new_file;
+        items="${items} ${files} ${new_file}";;
+  esac
+  shift
+done
+
+readchar=1
+while [ $readchar -ne 0 ]
+do
+  echo -n -e "\r" $items "- "
+  stty cbreak         # or stty raw
+  readchar=`dd if=/dev/tty bs=1 count=1 2>/dev/null`
+  stty -cbreak
+  if [ $readchar -ne 0 ] ; then
+    play -f s -r 8000 -s w ${file[$readchar]} $dsp > /dev/null
+  fi
+done
+echo
diff --git a/codec2/script/playraw.sh b/codec2/script/playraw.sh
index 0b6a3f6c..f6033e20 100755
--- a/codec2/script/playraw.sh
+++ b/codec2/script/playraw.sh
@@ -1,2 +1,6 @@
 #!/bin/sh
-play -f s -r 8000 -s w $1
+# Plays a raw file
+# usage:
+#   playraw file.raw
+#   playraw file.raw -d /deve/dsp1
+play -f s -r 8000 -s w $1 $2 $3
diff --git a/codec2/script/wav2raw.sh b/codec2/script/wav2raw.sh
index ab2290b8..39c0f1ae 100755
--- a/codec2/script/wav2raw.sh
+++ b/codec2/script/wav2raw.sh
@@ -1,2 +1,3 @@
 #!/bin/sh
+# Converts wave files to raw (headerless) files
 sox $1 -t raw $2
diff --git a/codec2/src/Makefile b/codec2/src/Makefile
index b83afdb0..164844b1 100644
--- a/codec2/src/Makefile
+++ b/codec2/src/Makefile
@@ -2,7 +2,8 @@ CC=gcc
 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
+SINEDEC_OBJ = sinedec.o globals.o initenc.o initdec.o four1.o synth.o \
+              quantise.o lpc.o dump.o
 
 all: sinenc sinedec
 
diff --git a/codec2/src/defines.h b/codec2/src/defines.h
index b524b80c..13898e2e 100644
--- a/codec2/src/defines.h
+++ b/codec2/src/defines.h
@@ -26,6 +26,9 @@
   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
+#ifndef __SINE__
+#define __SINE__
+
 /*---------------------------------------------------------------------------*\
                                                                              
 				DEFINES                                       
@@ -71,7 +74,6 @@ typedef struct {
   float imag;
 } COMP;
 
-
 /* Structure to hold unquantised model parameters for one frame */
 
 typedef struct {
@@ -81,3 +83,5 @@ typedef struct {
   float A[MAX_AMP];	/* average magnitude/unit frequency samples */
   float phi[MAX_AMP];	/* phase of each harmonic */
 } MODEL;
+
+#endif
diff --git a/codec2/src/dump.c b/codec2/src/dump.c
new file mode 100644
index 00000000..0deb15ca
--- /dev/null
+++ b/codec2/src/dump.c
@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+                                                                             
+  FILE........: dump.c
+  AUTHOR......: David Rowe                                                          
+  DATE CREATED: 25/8/09                                                       
+                                                                             
+  Routines to dump data to text files for Octave analysis.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License version 2, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#include "dump.h"
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+static int dumpon = 0;
+
+static FILE *fsn;
+static FILE *fmodel;
+static FILE *fqmodel;
+static FILE *fpw;
+
+void dump_on(char prefix[]) {
+    char s[MAX_STR];
+
+    dumpon = 1;
+
+    sprintf(s,"%s_sn.txt", prefix);
+    fsn = fopen(s, "wt");
+    assert(fsn != NULL);
+
+    sprintf(s,"%s_model.txt", prefix);
+    fmodel = fopen(s, "wt");
+    assert(fmodel != NULL);
+
+    sprintf(s,"%s_qmodel.txt", prefix);
+    fqmodel = fopen(s, "wt");
+    assert(fqmodel != NULL);
+
+    sprintf(s,"%s_pw.txt", prefix);
+    fpw = fopen(s, "wt");
+    assert(fpw != NULL);
+}
+
+void dump_off(){
+    fclose(fsn);
+    fclose(fmodel);
+    fclose(fqmodel);
+    fclose(fpw);
+}
+
+void dump_Sn(float Sn[]) {
+    int i;
+
+    if (!dumpon) return;
+
+    /* split across two lines to avoid max line length problems */
+    /* reconstruct in Octave */
+
+    for(i=0; i<AW_ENC/2; i++)
+	fprintf(fsn,"%f\t",Sn[i]);
+    fprintf(fsn,"\n");    
+    for(i=AW_ENC/2; i<AW_ENC; i++)
+	fprintf(fsn,"%f\t",Sn[i]);
+    fprintf(fsn,"\n");    
+}
+
+void dump_model(MODEL *model) {
+    int l;
+
+    if (!dumpon) return;
+
+    fprintf(fmodel,"%f\t%d\t", model->Wo, model->L);    
+    for(l=1; l<=model->L; l++)
+	fprintf(fmodel,"%f\t",model->A[l]);
+    for(l=model->L+1; l<MAX_AMP; l++)
+	fprintf(fmodel,"0.0\t");
+    fprintf(fmodel,"\n");    
+}
+
+void dump_quantised_model(MODEL *model) {
+    int l;
+
+    if (!dumpon) return;
+
+    fprintf(fqmodel,"%f\t%d\t", model->Wo, model->L);    
+    for(l=1; l<=model->L; l++)
+	fprintf(fqmodel,"%f\t",model->A[l]);
+    for(l=model->L+1; l<MAX_AMP; l++)
+	fprintf(fqmodel,"0.0\t");
+    fprintf(fqmodel,"\n");    
+}
+
+void dump_Pw(COMP Pw[]) {
+    int i;
+
+    if (!dumpon) return;
+
+    for(i=0; i<FFT_DEC/2; i++)
+	fprintf(fpw,"%f\t",Pw[i].real);
+    fprintf(fpw,"\n");    
+}
+
+
diff --git a/codec2/src/dump.h b/codec2/src/dump.h
new file mode 100644
index 00000000..61e41b2b
--- /dev/null
+++ b/codec2/src/dump.h
@@ -0,0 +1,39 @@
+/*---------------------------------------------------------------------------*\
+                                                                             
+  FILE........: dump.h
+  AUTHOR......: David Rowe                                                          
+  DATE CREATED: 25/8/09                                                       
+                                                                             
+  Routines to dump data to text files for Octave analysis.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License version 2, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#ifndef __DUMP__
+#define __DUMP__
+
+#include "sine.h"
+
+void dump_on(char filename_prefix[]);
+void dump_off();
+void dump_Sn(float Sn[]);
+void dump_model(MODEL *m);
+void dump_quantised_model(MODEL *m);
+void dump_Pw(COMP Pw[]);
+
+#endif
diff --git a/codec2/src/quantise.c b/codec2/src/quantise.c
new file mode 100644
index 00000000..e995ba2c
--- /dev/null
+++ b/codec2/src/quantise.c
@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+                                                                             
+  FILE........: quantise.c
+  AUTHOR......: David Rowe                                                          
+  DATE CREATED: 31/5/92                                                       
+                                                                             
+  Quantisation functions for the sinusoidal coder.  
+                                                                             
+\*---------------------------------------------------------------------------*/
+
+/*
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License version 2, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#include "sine.h"
+#include "quantise.h"
+#include "lpc.h"
+#include "dump.h"
+
+#define MAX_ORDER 20
+
+/*---------------------------------------------------------------------------*\
+									      
+  lpc_model_amplitudes
+
+  Derive a LPC model for amplitude samples then estimate amplitude samples
+  from this model with optional LSP quantisation.
+
+  Returns the spectral distortion for this frame.
+
+\*---------------------------------------------------------------------------*/
+
+float lpc_model_amplitudes(
+  float  Sn[],			/* Input frame of speech samples */
+  MODEL *model,			/* sinusoidal model parameters */
+  int    order,                 /* LPC model order */
+  int    lsp                    /* optional LSP quantisation if non-zero */
+)
+{
+  float Wn[AW_ENC];
+  float R[MAX_ORDER+1];
+  float ak[MAX_ORDER+1];
+  float E;
+  int   i;
+  float sd;			/* spectral distortion for this frame */
+
+  for(i=0; i<AW_ENC; i++)
+      Wn[i] = Sn[i]*w[i];
+  autocorrelate(Wn,R,AW_ENC,order);
+  levinson_durbin(R,ak,order);
+  E = 0.0;
+  for(i=0; i<=order; i++)
+      E += ak[i]*R[i];
+  
+  aks_to_M2(ak,order,model,E,&sd);   /* {ak} -> {Am} LPC decode */
+
+  return sd;
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                         
+   aks_to_M2()                                                             
+                                                                         
+   Transforms the linear prediction coefficients to spectral amplitude    
+   samples.  This function determines A(m) from the average energy per    
+   band using an FFT.                                                     
+                                                                        
+\*---------------------------------------------------------------------------*/
+
+void aks_to_M2(
+  float ak[],	/* LPC's */
+  int   order,
+  MODEL *model,	/* sinusoidal model parameters for this frame */
+  float E,	/* energy term */
+  float *sd	/* spectral distortion for this frame in dB */
+)
+{
+  COMP Pw[FFT_DEC];	/* power spectrum */
+  int i,m;		/* loop variables */
+  int am,bm;		/* limits of current band */
+  float r;		/* no. rads/bin */
+  float Em;		/* energy in band */
+  float Am;		/* spectral amplitude sample */
+  float noise;
+
+  r = TWO_PI/(FFT_DEC);
+
+  /* Determine DFT of A(exp(jw)) --------------------------------------------*/
+
+  for(i=0; i<FFT_DEC; i++) {
+    Pw[i].real = 0.0;
+    Pw[i].imag = 0.0;
+  }
+
+  for(i=0; i<=order; i++)
+    Pw[i].real = ak[i];
+  four1(&Pw[-1].imag,FFT_DEC,1);
+
+  /* Determine power spectrum P(w) = E/(A(exp(jw))^2 ------------------------*/
+
+  for(i=0; i<FFT_DEC/2; i++)
+    Pw[i].real = E/(Pw[i].real*Pw[i].real + Pw[i].imag*Pw[i].imag);
+  dump_Pw(Pw);
+
+  /* Determine magnitudes by linear interpolation of P(w) -------------------*/
+
+  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);
+    Em = 0.0;
+
+    for(i=am; i<bm; i++)
+      Em += Pw[i].real;
+    Am = sqrt(Em);
+
+    noise += pow(log10(Am/model->A[m]),2.0);
+    model->A[m] = Am;
+  }
+  *sd = 20.0*sqrt(noise/model->L);
+
+}
diff --git a/codec2/src/quantise.h b/codec2/src/quantise.h
new file mode 100644
index 00000000..00583ab9
--- /dev/null
+++ b/codec2/src/quantise.h
@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+                                                                             
+  FILE........: quantise.h
+  AUTHOR......: David Rowe                                                          
+  DATE CREATED: 31/5/92                                                       
+                                                                             
+  Quantisation functions for the sinusoidal coder.  
+                                                                             
+\*---------------------------------------------------------------------------*/
+
+/*
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License version 2, as
+  published by the Free Software Foundation.  This program is
+  distributed in the hope that it will be useful, but WITHOUT ANY
+  WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
+  License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program; if not, write to the Free Software
+  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+*/
+
+#ifndef __QUANTISE__
+#define __QUANTISE__
+
+#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);
+
+#endif
diff --git a/codec2/src/sinedec.c b/codec2/src/sinedec.c
index fece4e8a..ea0de97d 100644
--- a/codec2/src/sinedec.c
+++ b/codec2/src/sinedec.c
@@ -28,6 +28,8 @@
 
 #include <string.h>
 #include "sine.h"
+#include "quantise.h"
+#include "dump.h"
 
 /*---------------------------------------------------------------------------*\
                                                                              
@@ -68,11 +70,17 @@ int main(int argc, char *argv[])
   int i;		/* loop variable */
   int length;		/* number of frames so far */
 
-  char out_file[MAX_STR];
-  int arg;
+  char  out_file[MAX_STR];
+  int   arg;
+  float sd;
+  float sum_sd;
+
+  int lpc_model, order;
+  int dump;
 
   if (argc < 3) {
-    printf("usage: sinedec InputFile ModelFile [-o OutputFile]\n");
+    printf("usage: sinedec InputFile ModelFile [-o OutputFile] [-o lpc Order]\n");
+    printf("       [--dump DumpFilePrefix]\n");
     exit(0);
   }
 
@@ -116,6 +124,21 @@ int main(int argc, char *argv[])
   else
     length = 32000;
 
+  lpc_model = 0;
+  if ((arg = switch_present("--lpc",argc,argv))) {
+      lpc_model = 1;
+      order = atoi(argv[arg+1]);
+      if ((order < 4) || (order > 20)) {
+        printf("Error in lpc order: %d\n", order);
+        exit(1);
+      }	  
+  }
+
+  dump = switch_present("--dump",argc,argv);
+  if (dump) {
+      dump_on(argv[dump+1]);
+  }
+
   /* Initialise ------------------------------------------------------------*/
 
   init_decoder();
@@ -126,7 +149,8 @@ int main(int argc, char *argv[])
   /* Main loop ------------------------------------------------------------*/
 
   frames = 0;
-  while(fread(&model,sizeof(model),1,fmodel) /*&& frames < 1200*/) {
+  sum_sd = 0;
+  while(fread(&model,sizeof(model),1,fmodel)) {
     frames++;
 
     /* Read input speech */
@@ -136,6 +160,18 @@ int main(int argc, char *argv[])
       Sn[i] = Sn[i+N];
     for(i=0; i<N; i++)
       Sn[i+N+AW_ENC/2] = buf[i];
+    dump_Sn(Sn);
+ 
+    dump_model(&model);
+
+    /* optional LPC model amplitudes */
+
+    if (lpc_model) {
+	sd = lpc_model_amplitudes(Sn, &model, order, 0);
+	sum_sd += sd;
+    }
+
+    dump_quantised_model(&model);
 
     /* Synthesise speech */
 
@@ -154,13 +190,18 @@ int main(int argc, char *argv[])
 		buf[i] = Sn_[i];
 	}
 	fwrite(buf,sizeof(short),N,fout);
-    }
-
+    }    
   }
 
   if (fout != NULL)
     fclose(fout);
 
+  if (lpc_model)
+      printf("sd av = %5.2f dB\n", sum_sd/frames);
+
+  if (dump)
+      dump_off();
+
   return 0;
 }