--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ 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);
+}
+
--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ 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
--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ 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;
+}
projects / freetel-svn-tracking.git / commitdiff