From 1d3b30b3f9c4aa1a249ed867f96a91490f30b892 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Mon, 24 Aug 2009 06:21:25 +0000
Subject: [PATCH] added genres and some support octave and shell scripts

git-svn-id: https://svn.code.sf.net/p/freetel/code@29 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2/octave/.gitignore |   0
 codec2/octave/pl.m       |  42 +++++++
 codec2/octave/pl2.m      |  50 ++++++++
 codec2/script/playraw.sh |   2 +
 codec2/script/wav2raw.sh |   2 +
 codec2/src/lpc.c         | 250 +++++++++++++++++++++++++++++++++++++++
 codec2/src/lpc.h         |  40 +++++++
 codec2/unittest/Makefile |   9 ++
 codec2/unittest/genres.c |  92 ++++++++++++++
 9 files changed, 487 insertions(+)
 delete mode 100644 codec2/octave/.gitignore
 create mode 100644 codec2/octave/pl.m
 create mode 100644 codec2/octave/pl2.m
 create mode 100755 codec2/script/playraw.sh
 create mode 100755 codec2/script/wav2raw.sh
 create mode 100644 codec2/src/lpc.c
 create mode 100644 codec2/src/lpc.h
 create mode 100644 codec2/unittest/Makefile
 create mode 100644 codec2/unittest/genres.c

diff --git a/codec2/octave/.gitignore b/codec2/octave/.gitignore
deleted file mode 100644
index e69de29b..00000000
diff --git a/codec2/octave/pl.m b/codec2/octave/pl.m
new file mode 100644
index 00000000..dc33beb4
--- /dev/null
+++ b/codec2/octave/pl.m
@@ -0,0 +1,42 @@
+% Copyright David Rowe 2009
+% This program is distributed under the terms of the GNU General Public License 
+% Version 2
+
+function pl(samname1, start_sam, end_sam, pngname)
+  
+  fs=fopen(samname1,"rb");
+  s=fread(fs,Inf,"short");
+
+  st = 1;
+  en = length(s);
+  if (nargin >= 2)
+    st = start_sam;
+  endif
+  if (nargin >= 3)
+    en = end_sam;
+  endif
+
+  figure(1);
+  clg;
+  plot(s(st:en));
+  axis([1 en-st min(s) max(s)]);
+ 
+  if (nargin == 4)
+
+    % small image
+
+    __gnuplot_set__ terminal png size 420,300
+    ss = sprintf("__gnuplot_set__ output \"%s.png\"", pngname);
+    eval(ss)
+    replot;
+
+    % larger image
+
+    __gnuplot_set__ terminal png size 800,600
+    ss = sprintf("__gnuplot_set__ output \"%s_large.png\"", pngname);
+    eval(ss)
+    replot;
+
+  endif
+
+endfunction
diff --git a/codec2/octave/pl2.m b/codec2/octave/pl2.m
new file mode 100644
index 00000000..e18ca580
--- /dev/null
+++ b/codec2/octave/pl2.m
@@ -0,0 +1,50 @@
+% Copyright David Rowe 2009
+% This program is distributed under the terms of the GNU General Public License 
+% Version 2
+
+function pl2(samname1, samname2, start_sam, end_sam, pngname)
+  
+  fs1=fopen(samname1,"rb");
+  s1=fread(fs1,Inf,"short");
+  fs2=fopen(samname2,"rb");
+  s2=fread(fs2,Inf,"short");
+
+  st = 1;
+  en = length(s1);
+  if (nargin >= 3)
+    st = start_sam;
+  endif
+  if (nargin >= 4)
+    en = end_sam;
+  endif
+
+  figure(1);
+  clg;
+  subplot(211);
+  l1 = strcat("r;",samname1,";");
+  plot(s1(st:en), l1);
+  axis([1 en-st min(s1(st:en)) max(s1(st:en))]);
+  subplot(212);
+  l2 = strcat("r;",samname2,";");
+  plot(s2(st:en),l2);
+  axis([1 en-st min(s1(st:en)) max(s1(st:en))]);
+ 
+  if (nargin == 5)
+
+    % small image
+
+    __gnuplot_set__ terminal png size 420,300
+    s = sprintf("__gnuplot_set__ output \"%s.png\"", pngname);
+    eval(s)
+    replot;
+
+    % larger image
+
+    __gnuplot_set__ terminal png size 800,600
+    s = sprintf("__gnuplot_set__ output \"%s_large.png\"", pngname);
+    eval(s)
+    replot;
+
+  endif
+
+endfunction
diff --git a/codec2/script/playraw.sh b/codec2/script/playraw.sh
new file mode 100755
index 00000000..0b6a3f6c
--- /dev/null
+++ b/codec2/script/playraw.sh
@@ -0,0 +1,2 @@
+#!/bin/sh
+play -f s -r 8000 -s w $1
diff --git a/codec2/script/wav2raw.sh b/codec2/script/wav2raw.sh
new file mode 100755
index 00000000..ab2290b8
--- /dev/null
+++ b/codec2/script/wav2raw.sh
@@ -0,0 +1,2 @@
+#!/bin/sh
+sox $1 -t raw $2
diff --git a/codec2/src/lpc.c b/codec2/src/lpc.c
new file mode 100644
index 00000000..62bbed17
--- /dev/null
+++ b/codec2/src/lpc.c
@@ -0,0 +1,250 @@
+/*---------------------------------------------------------------------------*\
+                                                                           
+  FILE........: lpc.c                                                              
+  AUTHOR......: David Rowe                                                      
+  DATE CREATED: 30/9/90                                                   
+                                                                          
+  Linear Prediction functions written in C.                                
+                                                                          
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 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.
+*/
+
+#define LPC_MAX_N 512		/* maximum no. of samples in frame */
+#define MAX_ORDER 20		/* maximum LPC order */
+#define PI 3.141592654		/* mathematical constant */
+
+#include <lpc.h>
+#include <math.h>
+
+/*---------------------------------------------------------------------------*\
+                                                                         
+  hanning_window()                                                        
+                                                                         
+  Hanning windows a frame of speech samples.                              
+                                                                           
+\*---------------------------------------------------------------------------*/
+
+void hanning_window(
+  float Sn[],	/* input frame of speech samples */
+  float Wn[],	/* output frame of windowed samples */
+  int Nsam	/* number of samples */
+)
+{
+  int i;	/* loop variable */
+
+  for(i=0; i<Nsam; i++)
+    Wn[i] = Sn[i]*(0.5 - 0.5*cos(2*PI*(float)i/(Nsam-1)));
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                           
+  autocorrelate()                                                          
+                                                                          
+  Finds the first P autocorrelation values of an array of windowed speech 
+  samples Sn[].                                                            
+                                                                          
+\*---------------------------------------------------------------------------*/
+
+void autocorrelate(
+  float Sn[],	/* frame of Nsam windowed speech samples */
+  float Rn[],	/* array of P+1 autocorrelation coefficients */
+  int Nsam,	/* number of windowed samples to use */
+  int order	/* order of LPC analysis */
+)
+{
+  int i,j;	/* loop variables */
+
+  for(j=0; j<order+1; j++) {
+    Rn[j] = 0.0;
+    for(i=0; i<Nsam-j; i++)
+      Rn[j] += Sn[i]*Sn[i+j];
+  }
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                            
+  levinson_durbin()                                                        
+                                                                           
+  Given P+1 autocorrelation coefficients, finds P Linear Prediction Coeff. 
+  (LPCs) where P is the order of the LPC all-pole model. The Levinson-Durbin
+  algorithm is used, and is described in:                                   
+                                                                           
+    J. Makhoul                                                               
+    "Linear prediction, a tutorial review"                                   
+    Proceedings of the IEEE                                                
+    Vol-63, No. 4, April 1975                                               
+                                                                             
+\*---------------------------------------------------------------------------*/
+
+void levinson_durbin(
+  float R[],		/* order+1 autocorrelation coeff */
+  float lpcs[],		/* order+1 LPC's */
+  int order		/* order of the LPC analysis */
+)
+{
+  float E[MAX_ORDER+1];
+  float k[MAX_ORDER+1];
+  float a[MAX_ORDER+1][MAX_ORDER+1];
+  float sum;
+  int i,j;				/* loop variables */
+
+  E[0] = R[0];				/* Equation 38a, Makhoul */
+
+  for(i=1; i<=order; i++) {
+    sum = 0.0;
+    for(j=1; j<=i-1; j++)
+      sum += a[i-1][j]*R[i-j];
+    k[i] = -1.0*(R[i] + sum)/E[i-1];	/* Equation 38b, Makhoul */
+    if (fabs(k[i]) > 1.0)
+      k[i] = 0.0;
+
+    a[i][i] = k[i];
+
+    for(j=1; j<=i-1; j++)
+      a[i][j] = a[i-1][j] + k[i]*a[i-1][i-j];	/* Equation 38c, Makhoul */
+
+    E[i] = (1-k[i]*k[i])*E[i-1];		/* Equation 38d, Makhoul */
+  }
+
+  for(i=1; i<=order; i++)
+    lpcs[i] = a[order][i];
+  lpcs[0] = 1.0;  
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                             
+  inverse_filter()                                                          
+                                                                           
+  Inverse Filter, A(z).  Produces an array of residual samples from an array
+  of input samples and linear prediction coefficients.                      
+                                                                             
+  The filter memory is stored in the first order samples of the input array.
+                                                                             
+\*---------------------------------------------------------------------------*/
+
+void inverse_filter(
+  float Sn[],	/* Nsam input samples */
+  float a[],	/* LPCs for this frame of samples */
+  int Nsam,	/* number of samples */
+  float res[],	/* Nsam residual samples */
+  int order	/* order of LPC */
+)
+{
+  int i,j;	/* loop variables */
+
+  for(i=0; i<Nsam; i++) {
+    res[i] = 0.0;
+    for(j=0; j<=order; j++)
+      res[i] += Sn[i-j]*a[j];
+  }    
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                            
+ synthesis_filter()                                                        
+                                                                           
+ C version of the Speech Synthesis Filter, 1/A(z).  Given an array of   
+ residual or excitation samples, and the the LP filter coefficients, this
+ function will produce an array of speech samples.  This filter structure is
+ IIR.                                                                       
+                                                                           
+ The synthesis filter has memory as well, this is treated in the same way 
+ as the memory for the inverse filter (see inverse_filter() notes above). 
+ The difference is that the memory for the synthesis filter is stored in  
+ the output array, wheras the memory of the inverse filter is stored in the
+ input array.                                                              
+                                                                           
+ Note: the calling function must update the filter memory.                
+                                                                             
+\*---------------------------------------------------------------------------*/
+
+void synthesis_filter(
+  float res[],	/* Nsam input residual (excitation) samples */
+  float a[],	/* LPCs for this frame of speech samples */
+  int Nsam,	/* number of speech samples */
+  int order,	/* LPC order */
+  float Sn_[]	/* Nsam output synthesised speech samples */
+)
+{
+  int i,j;	/* loop variables */
+
+  /* Filter Nsam samples */
+
+  for(i=0; i<Nsam; i++) {
+    Sn_[i] = res[i]*a[0];
+    for(j=1; j<=order; j++)
+      Sn_[i] -= Sn_[i-j]*a[j];
+  }
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                            
+  find_aks()                                                                 
+                                                                            
+  This function takes a frame of samples, and determines the linear           
+  prediction coefficients for that frame of samples.                         
+                                                                            
+\*---------------------------------------------------------------------------*/
+
+void find_aks(
+  float Sn[],	/* Nsam samples with order sample memory */
+  float a[],	/* order+1 LPCs with first coeff 1.0 */
+  int Nsam,	/* number of input speech samples */
+  int order,	/* order of the LPC analysis */
+  float *E	/* residual energy */
+)
+{
+  float Wn[LPC_MAX_N];	/* windowed frame of Nsam speech samples */
+  float R[MAX_ORDER+1];	/* order+1 autocorrelation values of Sn[] */
+  int i;
+
+  hanning_window(Sn,Wn,Nsam);
+
+  autocorrelate(Wn,R,Nsam,order);
+  levinson_durbin(R,a,order);
+
+  *E = 0.0;
+  for(i=0; i<=order; i++)
+    *E += a[i]*R[i];
+  if (*E < 0.0)
+    *E = 1E-12;
+}
+
+/*---------------------------------------------------------------------------*\
+                                                                            
+  weight()                                                                  
+                                                                          
+  Weights a vector of LPCs.						   
+                                                                          
+\*---------------------------------------------------------------------------*/
+
+void weight(
+  float ak[],	/* vector of order+1 LPCs */
+  float gamma,	/* weighting factor */
+  int order,	/* num LPCs (excluding leading 1.0) */
+  float akw[]	/* weighted vector of order+1 LPCs */
+)
+{
+  int i;
+  
+  for(i=0; i<=order; i++)
+    akw[i] = ak[i]*pow(gamma,(float)i);
+}
+    
diff --git a/codec2/src/lpc.h b/codec2/src/lpc.h
new file mode 100644
index 00000000..dd133340
--- /dev/null
+++ b/codec2/src/lpc.h
@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+                                                                           
+  FILE........: lpc.h                                                            
+  AUTHOR......: David Rowe                                                      
+  DATE CREATED: 24/8/09                                                   
+                                                                          
+  Linear Prediction functions written in C.                                
+                                                                          
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 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.
+*/
+
+#ifndef __LPC__
+#define __LPC__
+
+void hanning_window(float Sn[],	float Wn[], int Nsam);
+void autocorrelate(float Sn[], float Rn[], int Nsam, int order);
+void levinson_durbin(float R[],	float lpcs[], int order);
+void inverse_filter(float Sn[], float a[], int Nsam, float res[], int order);
+void synthesis_filter(float res[], float a[], int Nsam,	int order, float Sn_[]);
+void find_aks(float Sn[], float a[], int Nsam, int order, float *E);
+void weight(float ak[],	float gamma, int order,	float akw[]);
+
+#endif
diff --git a/codec2/unittest/Makefile b/codec2/unittest/Makefile
new file mode 100644
index 00000000..898c13cc
--- /dev/null
+++ b/codec2/unittest/Makefile
@@ -0,0 +1,9 @@
+CFLAGS = -I../src -Wall -g
+
+all: genres
+
+genres: genres.o ../src/lpc.o
+	gcc $(CFLAGS) -o genres genres.o ../src/lpc.o -lm
+
+%.o : %.c
+	$(CC) -c $(CFLAGS) $< -o $@
diff --git a/codec2/unittest/genres.c b/codec2/unittest/genres.c
new file mode 100644
index 00000000..73770253
--- /dev/null
+++ b/codec2/unittest/genres.c
@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+                                                                           
+  FILE........: genres.c   
+  AUTHOR......: David Rowe                                                      
+  DATE CREATED: 24/8/09                                                   
+                                                                          
+  Generates a file of LPC residual samples from original speech.
+                                                                          
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2009 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.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <lpc.h>
+
+#define N 160
+#define P 10
+
+int main(int argc, char *argv[])
+{
+  FILE *fin,*fres;      /* input and output files */
+  short buf[N];         /* buffer of 16 bit speech samples */
+  float Sn[P+N];        /* input speech samples */
+  float res[N];         /* residual after LPC filtering */
+  float E;
+  float ak[P+1];        /* LP coeffs */
+
+  int frames;           /* frames processed so far */
+  int i;                /* loop variables */
+
+  if (argc < 3) {
+    printf("usage: %s InputFile ResidualFile\n", argv[0]);
+    exit(0);
+  }
+
+  /* Open files */
+
+  if ((fin = fopen(argv[1],"rb")) == NULL) {
+    printf("Error opening input file: %s\n",argv[1]);
+    exit(0);
+  }
+
+  if ((fres = fopen(argv[2],"wb")) == NULL) {
+    printf("Error opening output residual file: %s\n",argv[2]);
+    exit(0);
+  }
+
+  /* Initialise */
+
+  frames = 0;
+  for(i=0; i<P; i++) {
+    Sn[i] = 0.0;
+  }
+
+  /* Main loop */
+
+  while( (fread(buf,sizeof(short),N,fin)) == N) {
+    frames++;
+    for(i=0; i<N; i++)
+      Sn[P+i] = (float)buf[i];
+
+    /* Determine {ak} and filter to find residual */
+
+    find_aks(&Sn[P], ak, N, P, &E);
+    inverse_filter(&Sn[P], ak, N, res, P);
+    for(i=0; i<N; i++)
+      buf[i] = (short)res[i];
+    fwrite(buf,sizeof(short),N,fres);
+  }
+
+  fclose(fin);
+  fclose(fres);
+
+  return 0;
+}
-- 
2.25.1