CC=gcc
CFLAGS=-g -Wall -I. -I../src -Wall -g -DFLOATING_POINT -DVAR_ARRAYS
-SINENC_OBJ = sinenc.o globals.o initenc.o four1.o refine.o spec.o dump.o
-SINEDEC_OBJ = sinedec.o globals.o initenc.o initdec.o four1.o synth.o \
- quantise.o lpc.o dump.o refine.o lsp.o phase.o postfilter.o \
- interp.o nlp.o
+SINENC_OBJ = sine.o four1.o sinenc.o dump.o
+SINEDEC_OBJ = sine.o four1.o dump.o quantise.o lpc.o lsp.o phase.o \
+ postfilter.o interp.o sinedec.o
all: sinenc sinedec
../unittest/tnlp ../raw/$1.raw ../unittest/$1_nlp.p
../src/sinenc ../raw/$1.raw $1.mdl 300 ../unittest/$1_nlp.p
../src/sinedec ../raw/$1.raw $1.mdl -o $1_uq.raw
-../src/sinedec ../raw/$1.raw $1.mdl --phase 0 -o $1_phase0.raw --postfilter
+../src/sinedec ../raw/$1.raw $1.mdl --phase0 -o $1_phase0.raw --postfilter
../src/sinedec ../raw/$1.raw $1.mdl --lpc 10 -o $1_lpc10.raw
../src/sinedec ../raw/$1.raw $1.mdl --lpc 10 --lsp -o $1_lsp.raw
-../src/sinedec ../raw/$1.raw $1.mdl --phase 0 --lpc 10 -o $1_phase0_lpc10.raw --postfilter
-../src/sinedec ../raw/$1.raw $1.mdl --phase 0 --lpc 10 --lsp -o $1_phase0_lsp.raw --postfilter
+../src/sinedec ../raw/$1.raw $1.mdl --phase0 --lpc 10 -o $1_phase0_lpc10.raw --postfilter
+../src/sinedec ../raw/$1.raw $1.mdl --phase0 --lpc 10 --lsp -o $1_phase0_lsp.raw --postfilter
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
-#ifndef __SINE__
-#define __SINE__
+#ifndef __DEFINES__
+#define __DEFINES__
/*---------------------------------------------------------------------------*\
/* General defines */
-#define N 80 /* number of samples per frame */
-#define MAX_AMP 80 /* maximum number of harmonics */
-#define PI 3.141592654 /* mathematical constant */
-#define TWO_PI 6.283185307 /* mathematical constant */
-#define FC 4000 /* cut-off frequency */
+#define N 80 /* number of samples per frame */
+#define MAX_AMP 80 /* maximum number of harmonics */
+#define PI 3.141592654 /* mathematical constant */
+#define TWO_PI 6.283185307 /* mathematical constant */
+#define FS 8000 /* sample rate in Hz */
+#define MAX_STR 256 /* maximum string size */
+
+#define NW 279 /* analysis window size */
+#define FFT_ENC 512 /* size of FFT used for encoder */
+#define FFT_DEC 512 /* size of FFT used in decoder */
+#define TW 40 /* Trapezoidal synthesis window overlap */
+#define V_THRESH 4.0 /* voicing threshold in dB */
+#define LPC_MAX 20 /* maximum LPC order */
+#define PHASE_LPC 10 /* maximum LPC order */
/* Pitch estimation defines */
-#define M 320 /* pitch analysis frame size */
-#define P_MIN 20 /* minimum pitch */
-#define P_MAX 160 /* maximum pitch */
-
-/* Encoder defines */
-
-#define NW 279 /* analysis window size */
-#define FFT_ENC 512 /* size of FFT used for encoder analysis */
-
-/* Decoder defines */
-
-#define AW_DEC 160 /* number of samples in synthesis window */
-#define FFT_DEC 512 /* number of points in DFT */
-#define TW 40 /* Trapezoidal synthesis window overlap */
-#define MAX_STR 256
+#define M 320 /* pitch analysis frame size */
+#define P_MIN 20 /* minimum pitch */
+#define P_MAX 160 /* maximum pitch */
/*---------------------------------------------------------------------------*\
float imag;
} COMP;
-/* Structure to hold unquantised model parameters for one frame */
+/* Structure to hold model parameters for one frame */
typedef struct {
- float Wo; /* fundamental frequency estimate in radians */
- int L; /* number of harmonics over the current frame */
- float v[MAX_AMP]; /* voicing measures */
- float A[MAX_AMP]; /* average magnitude/unit frequency samples */
- float phi[MAX_AMP]; /* phase of each harmonic */
+ float Wo; /* fundamental frequency estimate in radians */
+ int L; /* number of harmonics */
+ float v[MAX_AMP]; /* voicing measures (unused at present) */
+ float A[MAX_AMP]; /* amplitiude of each harmonic */
+ float phi[MAX_AMP]; /* phase of each harmonic */
} MODEL;
#endif
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
+#include "defines.h"
#include "dump.h"
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
+#include <math.h>
static int dumpon = 0;
fprintf(fqmodel,"\n");
}
-void dump_phase(float phase[]) {
+void dump_phase(float phase[], int L) {
int l;
char s[MAX_STR];
assert(fphase != NULL);
}
- for(l=1; l<=model.L; l++)
+ for(l=1; l<=L; l++)
fprintf(fphase,"%f\t",phase[l]);
- for(l=model.L+1; l<MAX_AMP; l++)
+ for(l=L+1; l<MAX_AMP; l++)
fprintf(fphase,"%f\t",0.0);
fprintf(fphase,"\n");
}
-void dump_phase_(float phase_[]) {
+void dump_phase_(float phase_[], int L) {
int l;
char s[MAX_STR];
assert(fphase_ != NULL);
}
- for(l=1; l<=model.L; l++)
+ for(l=1; l<=L; l++)
fprintf(fphase_,"%f\t",phase_[l]);
- for(l=model.L+1; l<MAX_AMP; l++)
+ for(l=L+1; l<MAX_AMP; l++)
fprintf(fphase_,"%f\t",0.0);
fprintf(fphase_,"\n");
}
#ifndef __DUMP__
#define __DUMP__
-#include "sine.h"
-
void dump_on(char filename_prefix[]);
void dump_off();
/* phase modelling */
void dump_snr(float snr);
-void dump_phase(float phase[]);
-void dump_phase_(float phase[]);
+void dump_phase(float phase[], int L);
+void dump_phase_(float phase[], int L);
/* NLP states */
+++ /dev/null
-/*---------------------------------------------------------------------------*\
-
- FILE: initenc.c
- AUTHOR: David Rowe
- DATE CREATED: 11/5/94
-
- Initialises sinusoidal speech encoder.
-
-\*---------------------------------------------------------------------------*/
-
-/*
- 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 "sine.h" /* sinusoidal header file */
-
-void init_encoder()
-{
- int i;
-
- frames = 0;
-
- /* Initialise sample buffer memories, 1 1 stop divide by zero errors
- and zero energy frames at the start */
-
- for(i=0; i<M; i++)
- Sn[i] = 1.0;
-
-}
-
-float make_window()
-{
- float m;
- COMP temp;
- int i,j;
-
- /*
- Generate Hamming window centered on M-sample pitch analysis window
-
- 0 M/2 M-1
- |-------------|-------------|
- |-------|-------|
- NW samples
-
- All our analysis/synthsis is centred on the M/2 sample.
- */
-
- m = 0.0;
- for(i=0; i<M/2-NW/2; i++)
- w[i] = 0.0;
- for(i=M/2-NW/2,j=0; i<M/2+NW/2; i++,j++) {
- w[i] = 0.5 - 0.5*cos(TWO_PI*j/(NW-1));
- m += w[i]*w[i];
- }
- for(i=M/2+NW/2; i<M; i++)
- w[i] = 0.0;
-
- /* Normalise - makes freq domain amplitude estimation straight
- forward */
-
- m = 1.0/sqrt(m*FFT_ENC);
- for(i=0; i<M; i++) {
- w[i] *= m;
- }
-
- /*
- Generate DFT of analysis window, used for later processing. Note
- we modulo FFT_ENC shift the time domain window w[], this makes the
- imaginary part of the DFT W[] equal to zero as the shifted w[] is
- even about the n=0 time axis if NW is odd. Having the imag part
- of the DFT W[] makes computation easier.
-
- 0 FFT_ENC-1
- |-------------------------|
-
- ----\ /----
- \ /
- \ / <- shifted version of window w[n]
- \ /
- \ /
- -------
-
- |---------| |---------|
- NW/2 NW/2
- */
-
- for(i=0; i<FFT_ENC; i++) {
- W[i].real = 0.0;
- W[i].imag = 0.0;
- }
- for(i=0; i<NW/2; i++)
- W[i].real = w[i+M/2];
- for(i=FFT_ENC-NW/2,j=M/2-NW/2; i<FFT_ENC; i++,j++)
- W[i].real = w[j];
-
- four1(&W[-1].imag,FFT_ENC,-1); /* "Numerical Recipes in C" FFT */
-
- /*
- Re-arrange W[] to be symmetrical about FFT_ENC/2. Makes later
- analysis convenient.
-
- Before:
-
-
- 0 FFT_ENC-1
- |----------|---------|
- __ _
- \ /
- \_______________/
-
- After:
-
- 0 FFT_ENC-1
- |----------|---------|
- ___
- / \
- ________/ \_______
-
- */
-
-
- for(i=0; i<FFT_ENC/2; i++) {
- temp.real = W[i].real;
- temp.imag = W[i].imag;
- W[i].real = W[i+FFT_ENC/2].real;
- W[i].imag = W[i+FFT_ENC/2].imag;
- W[i+FFT_ENC/2].real = temp.real;
- W[i+FFT_ENC/2].imag = temp.imag;
- }
-
- return(m);
-}
-
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
-#include "interp.h"
+#include <math.h>
#include <string.h>
+#include "defines.h"
+#include "interp.h"
+
/*---------------------------------------------------------------------------*\
interp()
phase to time-shift the model forward or backward N
samples. */
- memcpy(a, prev, sizeof(model));
- memcpy(b, next, sizeof(model));
+ memcpy(a, prev, sizeof(MODEL));
+ memcpy(b, next, sizeof(MODEL));
for(m=1; m<=a->L; m++) {
a->A[m] /= 2.0;
a->phi[m] += a->Wo*m*N;
#ifndef __INTERP__
#define __INTERP__
-#include "sine.h"
-
void interp(MODEL *prev, MODEL *next, MODEL *synth, MODEL *a, MODEL *b,
int *transition);
*/
#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 <assert.h>
#include <math.h>
+#include "defines.h"
+#include "lpc.h"
/*---------------------------------------------------------------------------*\
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 E[LPC_MAX+1];
+ float k[LPC_MAX+1];
+ float a[LPC_MAX+1][LPC_MAX+1];
float sum;
int i,j; /* loop variables */
)
{
float Wn[LPC_MAX_N]; /* windowed frame of Nsam speech samples */
- float R[MAX_ORDER+1]; /* order+1 autocorrelation values of Sn[] */
+ float R[LPC_MAX+1]; /* order+1 autocorrelation values of Sn[] */
int i;
- hanning_window(Sn,Wn,Nsam);
+ assert(order < LPC_MAX);
+ assert(Nsam < LPC_MAX_N);
+ hanning_window(Sn,Wn,Nsam);
autocorrelate(Wn,R,Nsam,order);
levinson_durbin(R,a,order);
}
}
- #ifdef POST_PROCESS_MBE
- best_f0 = post_process_mbe(Fw, pmin, pmax, gmax);
- #else
best_f0 = post_process_sub_multiples(Fw, pmin, pmax, gmax, gmax_bin, prev_Wo);
- #endif
/* Shift samples in buffer to make room for new samples */
return best_f0;
}
-/*---------------------------------------------------------------------------*\
-
- post_process_mbe()
-
- Use the MBE pitch estimation algorithm to evaluate pitch candidates. This
- works OK but the accuracy at low F0 is affected by NW, the analysis window
- size used for the DFT of the input speech Sw[]. Also favours high F0 in
- the presence of background noise which causes periodic artifacts in the
- synthesised speech.
-
-\*---------------------------------------------------------------------------*/
-
-float post_process_mbe(COMP Fw[], int pmin, int pmax, float gmax)
-{
- float candidate_f0;
- float f0,best_f0; /* fundamental frequency */
- float e,e_min; /* MBE cost function */
- int i;
- float e_hz[F0_MAX];
- int bin;
- float f0_min, f0_max;
- float f0_start, f0_end;
- COMP Sw_[FFT_ENC];
-
- f0_min = (float)SAMPLE_RATE/pmax;
- f0_max = (float)SAMPLE_RATE/pmin;
-
- /* Now look for local maxima. Each local maxima is a candidate
- that we test using the MBE pitch estimation algotithm */
-
- for(i=0; i<F0_MAX; i++)
- e_hz[i] = -1;
- e_min = 1E32;
- best_f0 = 50;
- for(i=PE_FFT_SIZE*DEC/pmax; i<=PE_FFT_SIZE*DEC/pmin; i++) {
- if ((Fw[i].real > Fw[i-1].real) && (Fw[i].real > Fw[i+1].real)) {
-
- /* local maxima found, lets test if it's big enough */
-
- if (Fw[i].real > T*gmax) {
-
- /* OK, sample MBE cost function over +/- 10Hz range in 2.5Hz steps */
-
- candidate_f0 = (float)i*SAMPLE_RATE/(PE_FFT_SIZE*DEC);
- f0_start = candidate_f0-20;
- f0_end = candidate_f0+20;
- if (f0_start < f0_min) f0_start = f0_min;
- if (f0_end > f0_max) f0_end = f0_max;
-
- for(f0=f0_start; f0<=f0_end; f0+= 2.5) {
- e = test_candidate_mbe(Sw, f0, Sw_);
- bin = floor(f0); assert((bin > 0) && (bin < F0_MAX));
- e_hz[bin] = e;
- if (e < e_min) {
- e_min = e;
- best_f0 = f0;
- }
- }
-
- }
- }
- }
- dump_e(e_hz);
-
- return best_f0;
-}
-
-/*---------------------------------------------------------------------------*\
-
- test_candidate_mbe()
-
- Returns the error of the MBE cost function for the input f0.
-
- Note: I think a lot of the operations below can be simplified as
- W[].imag = 0 and has been normalised such that den always equals 1.
-
-\*---------------------------------------------------------------------------*/
-
-float test_candidate_mbe(
- COMP Sw[],
- float f0,
- COMP Sw_[] /* DFT of all voiced synthesised signal for f0 */
- /* useful for debugging/dump file */
-)
-{
- int i,l,al,bl,m; /* loop variables */
- COMP Am; /* amplitude sample for this band */
- int offset; /* centers Hw[] about current harmonic */
- float den; /* denominator of Am expression */
- float error; /* accumulated error between originl and synthesised */
- float Wo; /* current "test" fundamental freq. */
- int L;
-
- for(i=0; i<FFT_ENC; i++) {
- Sw_[i].real = 0.0;
- Sw_[i].imag = 0.0;
- }
-
- L = floor((SAMPLE_RATE/2.0)/f0);
- Wo = f0*(2*PI/SAMPLE_RATE);
-
- error = 0.0;
-
- /* Just test across the harmonics in the first 1000 Hz (L/4) */
-
- for(l=1; l<=L/4; l++) {
- Am.real = 0.0;
- Am.imag = 0.0;
- den = 0.0;
- al = ceil((l - 0.5)*Wo*FFT_ENC/TWO_PI);
- bl = ceil((l + 0.5)*Wo*FFT_ENC/TWO_PI);
-
- /* Estimate amplitude of harmonic assuming harmonic is totally voiced */
-
- for(m=al; m<bl; m++) {
- offset = FFT_ENC/2 + m - l*Wo*FFT_ENC/TWO_PI + 0.5;
- Am.real += Sw[m].real*W[offset].real + Sw[m].imag*W[offset].imag;
- Am.imag += Sw[m].imag*W[offset].real - Sw[m].real*W[offset].imag;
- den += W[offset].real*W[offset].real + W[offset].imag*W[offset].imag;
- }
-
- Am.real = Am.real/den;
- Am.imag = Am.imag/den;
-
- /* Determine error between estimated harmonic and original */
-
- for(m=al; m<bl; m++) {
- offset = FFT_ENC/2 + m - l*Wo*FFT_ENC/TWO_PI + 0.5;
- Sw_[m].real = Am.real*W[offset].real - Am.imag*W[offset].imag;
- Sw_[m].imag = Am.real*W[offset].imag + Am.imag*W[offset].real;
- error += (Sw[m].real - Sw_[m].real)*(Sw[m].real - Sw_[m].real);
- error += (Sw[m].imag - Sw_[m].imag)*(Sw[m].imag - Sw_[m].imag);
- }
- }
-
- return error;
-}
-
-
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
-#include "sine.h"
+#include "defines.h"
#include "phase.h"
-#include "lpc.h"
+#include "four1.h"
+
#include <assert.h>
+#include <math.h>
#include <string.h>
+#include <stdlib.h>
#define VTHRESH 4.0
\*---------------------------------------------------------------------------*/
-void aks_to_H(model,aks,G,H, order)
-MODEL *model; /* model parameters */
-float aks[]; /* LPC's */
-float G; /* energy term */
-COMP H[]; /* complex LPC spectral samples */
-int order;
+void aks_to_H(
+ MODEL *model, /* model parameters */
+ float aks[], /* LPC's */
+ float G, /* energy term */
+ COMP H[], /* complex LPC spectral samples */
+ int order
+)
{
COMP Pw[FFT_DEC]; /* power spectrum */
int i,m; /* loop variables */
}
}
-/*---------------------------------------------------------------------------*\
-
- phase_model_first_order()
-
- Models the current frames phase samples {phi[m]} as a LPC synthesis
- filter excited by an impulse at time i_min with a complex gain min_Am.
- The phase of the complex gain min_Am is a constant phase term, and the
- impulse position i_min creates a phase term that varies linearly
- with frequency. We are therefore modelling the excitation phase as
- a 1st order polynomial function of mWo:
-
- Ex[m] = -m*Wo*i_min + arg(min_Am);
-
- and the harmonic phase is given by:
-
- phi[m] = -m*Wo*i_min + arg(min_Am) + arg[H[m]];
-
- where H[m] is the LPC spectra sampled at the mWo harmonic frequencies.
-
- This actually works out to be a pretty good model for voiced
- speech. The fit of the model (snr) is therefore used as a voicing
- measure. If the snr is less than a threshold, the frame is
- declared unvoiced all all the phases are randomised.
-
- Reference:
- [1] http://www.itr.unisa.edu.au/~steven/thesis/dgr.pdf Chapter 6
-
- NOTE: min_Am is a dumb name for the complex gain as {A} or A[m] is
- commonly used for the spectral magnitudes. TODO: change name to match
- thesis term G(mWo) (which of course clashes with LPC gain! AAAAAAHHH!).
-
-\*---------------------------------------------------------------------------*/
-
-float phase_model_first_order(
- float aks[], /* LPC coeffs for this frame */
- COMP H[], /* LPC filter freq doamin samples */
- float *n_min, /* pulse position for min error */
- COMP *minAm, /* complex gain for min error */
- int *voiced
-)
-{
- float G; /* LPC gain */
- int m;
-
- float E,Emin; /* current and minimum error so far */
- int P; /* current pitch period */
- COMP A[MAX_AMP]; /* harmonic samples */
- COMP Ex[MAX_AMP]; /* excitation samples */
- COMP A_[MAX_AMP]; /* synthesised harmonic samples */
- COMP Am; /* complex gain */
- COMP Em; /* error for m-th band */
- float den; /* energy of synthesised */
- float snr; /* snr of each excitation source */
- int Lmax;
- float n;
-
- Lmax = model.L;
-
- /* Construct target vector */
-
- sig = 0.0;
- for(m=1; m<=Lmax; m++) {
- A[m].real = model.A[m]*cos(model.phi[m]);
- A[m].imag = model.A[m]*sin(model.phi[m]);
- sig += model.A[m]*model.A[m];
- }
-
- /* Sample LPC model at harmonics */
-
- //#define NO_LPC_PHASE
- #ifdef NO_LPC_PHASE
- /* useful for testing with Sn[] an impulse train */
- for(m=1; m<=PHASE_LPC_ORD; m++)
- aks[m] = 0;
- #endif
- G = 1.0;
- aks_to_H(&model,aks,G,H,PHASE_LPC_ORD);
-
- /* Now attempt to fit impulse, by trying positions 0..P-1 */
-
- Emin = 1E32;
- P = floor(TWO_PI/model.Wo + 0.5);
- for(n=0; n<P; n+=0.25) {
-
- /* determine complex gain */
-
- Am.real = 0.0;
- Am.imag = 0.0;
- den = 0.0;
- for(m=1; m<=Lmax; m++) {
- Ex[m].real = cos(model.Wo*m*n);
- Ex[m].imag = sin(-model.Wo*m*n);
- A_[m].real = H[m].real*Ex[m].real - H[m].imag*Ex[m].imag;
- A_[m].imag = H[m].imag*Ex[m].real + H[m].real*Ex[m].imag;
- Am.real += A[m].real*A_[m].real + A[m].imag*A_[m].imag;
- Am.imag += A[m].imag*A_[m].real - A[m].real*A_[m].imag;
- den += A_[m].real*A_[m].real + A_[m].imag*A_[m].imag;
- }
-
- Am.real /= den;
- Am.imag /= den;
-
- /* determine error */
-
- E = 0.0;
- for(m=1; m<=Lmax; m++) {
- float new_phi;
-
- Em.real = A_[m].real*Am.real - A_[m].imag*Am.imag;
- Em.imag = A_[m].imag*Am.real + A_[m].real*Am.imag;
-
- new_phi = atan2(Em.imag, Em.real+1E-12);
- E += pow(model.A[m]*(cos(model.phi[m])-cos(new_phi)),2.0);
- E += pow(model.A[m]*(sin(model.phi[m])-sin(new_phi)),2.0);
- }
-
- if (E < Emin) {
- Emin = E;
- *n_min = n;
- minAm->real = Am.real;
- minAm->imag = Am.imag;
- }
-
- }
-
- snr = 10.0*log10(sig/Emin);
- if (snr > VTHRESH)
- *voiced = 1;
- else
- *voiced = 0;
-
- return snr;
-}
/*---------------------------------------------------------------------------*\
\*---------------------------------------------------------------------------*/
void phase_synth_zero_order(
- float snr, /* SNR from first order model */
- COMP H[], /* LPC spectra samples */
- float *ex_phase,/* excitation phase of fundamental */
- int voiced
+ MODEL *model,
+ float aks[],
+ int voiced,
+ float *ex_phase /* excitation phase of fundamental */
)
{
- int m,i;
+ int m;
float new_phi;
COMP Ex[MAX_AMP]; /* excitation samples */
COMP A_[MAX_AMP]; /* synthesised harmonic samples */
- float maxA;
+ COMP H[MAX_AMP]; /* LPC freq domain samples */
+ float G;
float jitter;
+ G = 1.0;
+ aks_to_H(model,aks,G,H,PHASE_LPC);
+
/*
Update excitation fundamental phase track, this sets the position
of each pitch pulse during voiced speech. After much experiment
I found that using just this frame Wo improved quality for UV
sounds compared to interpolating two frames Wo like this:
- ex_phase[0] += (*prev_Wo+model.Wo)*N/2;
+ ex_phase[0] += (*prev_Wo+mode->Wo)*N/2;
*/
- ex_phase[0] += (model.Wo)*N;
+ ex_phase[0] += (model->Wo)*N;
ex_phase[0] -= TWO_PI*floor(ex_phase[0]/TWO_PI + 0.5);
- for(m=1; m<=model.L; m++) {
+ for(m=1; m<=model->L; m++) {
/* generate excitation */
of each harmonic over at +/- 0.25 of a sample.
*/
jitter = 0.25*(1.0 - 2.0*rand()/RAND_MAX);
- Ex[m].real = cos(ex_phase[0]*m - jitter*model.Wo*m);
- Ex[m].imag = sin(ex_phase[0]*m - jitter*model.Wo*m);
+ Ex[m].real = cos(ex_phase[0]*m - jitter*model->Wo*m);
+ Ex[m].imag = sin(ex_phase[0]*m - jitter*model->Wo*m);
}
else {
/* modify sinusoidal phase */
new_phi = atan2(A_[m].imag, A_[m].real+1E-12);
- model.phi[m] = new_phi;
+ model->phi[m] = new_phi;
}
}
-
-/*---------------------------------------------------------------------------*\
-
- phase_synth_first_order()
-
- Synthesises phases based on SNR and the first oreder phase model
- parameters.
-
-\*---------------------------------------------------------------------------*/
-
-void phase_synth_first_order(
- float snr, /* SNR from first order model */
- COMP H[], /* LPC spectra samples */
- float n_min, /* best pulse position */
- COMP minAm, /* best complex gain */
- int voiced
-)
-{
- int Lrand;
- int m;
- float new_phi;
- COMP Ex[MAX_AMP]; /* excitation samples */
- COMP A_[MAX_AMP]; /* synthesised harmonic samples */
- COMP Tm;
-
- /* now modify sinusoidal model phase using phase model */
-
- for(m=1; m<=model.L; m++) {
-
- /* generate excitation */
-
- if (voiced) {
- Ex[m].real = cos(model.Wo*m*n_min);
- Ex[m].imag = sin(-model.Wo*m*n_min);
- }
- else {
- float phi = TWO_PI*(float)rand()/RAND_MAX;
- Ex[m].real = cos(phi);
- Ex[m].imag = sin(phi);
- }
-
- /* filter using LPC filter */
-
- A_[m].real = H[m].real*Ex[m].real - H[m].imag*Ex[m].imag;
- A_[m].imag = H[m].imag*Ex[m].real + H[m].real*Ex[m].imag;
-
- /* scale by complex gain (could have done this earlier at Ex[]
- stage) */
-
- Tm.real = A_[m].real*minAm.real - A_[m].imag*minAm.imag;
- Tm.imag = A_[m].imag*minAm.real + A_[m].real*minAm.imag;
-
- /* modify sinusoidal phase */
-
- new_phi = atan2(Tm.imag, Tm.real+1E-12);
- model.phi[m] = new_phi;
- }
-}
-
#ifndef __PHASE__
#define __PHASE__
-#define PHASE_LPC_ORD 10
-
-void aks_to_H(MODEL *model, float aks[], float G, COMP H[], int order);
-float phase_model_first_order(float aks[], COMP H[], float *n_min, COMP *min_Am,
- int *voiced);
-void phase_synth_zero_order(float snr, COMP H[], float *ex_phase, int voiced);
-void phase_synth_first_order(float snr, COMP H[], float n_min, COMP min_Am,
- int voiced);
+void phase_synth_zero_order(MODEL *model, float aks[], int voiced, float *ex_phase);
#endif
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
-#include "postfilter.h"
+
+#include "defines.h"
#include "dump.h"
+#include "postfilter.h"
/*---------------------------------------------------------------------------*\
This idea is rather experimental. Some potential problems that may
happen:
- 1/ If someone says "aaaaaaaahhhhhhhhh" willl background estimator track
+ 1/ If someone says "aaaaaaaahhhhhhhhh" will background estimator track
up to speech level? This would be a bad thing.
2/ If background noise suddenly dissapears from the source speech does
#ifndef __POSTFILTER__
#define __POSTFILTER__
-#include "sine.h"
-
void postfilter(MODEL *model, int voiced, float *bg_est);
#endif
#include <assert.h>
#include <ctype.h>
-#include "sine.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "defines.h"
+#include "dump.h"
#include "quantise.h"
#include "lpc.h"
#include "lsp.h"
-#include "dump.h"
+#include "four1.h"
-#define MAX_ORDER 20
#define LSP_DELTA1 0.01 /* grid spacing for LSP root searches */
#define MAX_CB 20 /* max number of codebooks */
)
{
int i;
- float dlsp[MAX_ORDER];
- float dlsp_[MAX_ORDER];
+ float dlsp[LPC_MAX];
+ float dlsp_[LPC_MAX];
dlsp[0] = lsp[0];
for(i=1; i<order; i++)
float lpc_model_amplitudes(
float Sn[], /* Input frame of speech samples */
+ float w[],
MODEL *model, /* sinusoidal model parameters */
int order, /* LPC model order */
int lsp_quant, /* optional LSP quantisation if non-zero */
)
{
float Wn[M];
- float R[MAX_ORDER+1];
+ float R[LPC_MAX+1];
float E;
int i,j;
float snr;
- float lsp[MAX_ORDER];
- float lsp_hz[MAX_ORDER];
- float lsp_[MAX_ORDER];
- int roots; /* number of LSP roots found */
+ float lsp[LPC_MAX];
+ float lsp_hz[LPC_MAX];
+ float lsp_[LPC_MAX];
+ int roots; /* number of LSP roots found */
int index;
float se;
int k,m;
float *cb;
- float wt[MAX_ORDER];
+ float wt[LPC_MAX];
for(i=0; i<M; i++)
Wn[i] = Sn[i]*w[i];
#ifndef __QUANTISE__
#define __QUANTISE__
-#include "sine.h"
-
void quantise_init();
-float lpc_model_amplitudes(float Sn[], MODEL *model, int order,int lsp,float ak[]);
+float lpc_model_amplitudes(float Sn[], float w[], MODEL *model, int order,int lsp,float ak[]);
void aks_to_M2(float ak[], int order, MODEL *model, float E, float *snr);
float get_gmin(void);
DATE CREATED: 27/5/94
Finds the DFT of the current speech input speech frame.
-
- INPUT.......: global float Sn[] speech samples
- OUTPUT......: global float Sw[] DFT of speech samples
\*---------------------------------------------------------------------------*/
-void dft_speech()
+void dft_speech(float Sn[], COMP Sw[])
{
int i;
--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ FILE........: sine.c
+ AUTHOR......: David Rowe
+ DATE CREATED: 119/8/2010
+
+ Sinusoidal analysis and synthesis functions.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ Copyright (C) 1990-2010 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.
+*/
+
+/*---------------------------------------------------------------------------*\
+
+ INCLUDES
+
+\*---------------------------------------------------------------------------*/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "defines.h"
+#include "sine.h"
+#include "four1.h"
+
+/*---------------------------------------------------------------------------*\
+
+ HEADERS
+
+\*---------------------------------------------------------------------------*/
+
+void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, float pstep);
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTIONS
+
+\*---------------------------------------------------------------------------*/
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: make_analysis_window
+ AUTHOR......: David Rowe
+ DATE CREATED: 11/5/94
+
+ Init function that generates the time domain analysis window and it's DFT.
+
+\*---------------------------------------------------------------------------*/
+
+void make_analysis_window(float w[],COMP W[])
+{
+ float m;
+ COMP temp;
+ int i,j;
+
+ /*
+ Generate Hamming window centered on M-sample pitch analysis window
+
+ 0 M/2 M-1
+ |-------------|-------------|
+ |-------|-------|
+ NW samples
+
+ All our analysis/synthsis is centred on the M/2 sample.
+ */
+
+ m = 0.0;
+ for(i=0; i<M/2-NW/2; i++)
+ w[i] = 0.0;
+ for(i=M/2-NW/2,j=0; i<M/2+NW/2; i++,j++) {
+ w[i] = 0.5 - 0.5*cos(TWO_PI*j/(NW-1));
+ m += w[i]*w[i];
+ }
+ for(i=M/2+NW/2; i<M; i++)
+ w[i] = 0.0;
+
+ /* Normalise - makes freq domain amplitude estimation straight
+ forward */
+
+ m = 1.0/sqrt(m*FFT_ENC);
+ for(i=0; i<M; i++) {
+ w[i] *= m;
+ }
+
+ /*
+ Generate DFT of analysis window, used for later processing. Note
+ we modulo FFT_ENC shift the time domain window w[], this makes the
+ imaginary part of the DFT W[] equal to zero as the shifted w[] is
+ even about the n=0 time axis if NW is odd. Having the imag part
+ of the DFT W[] makes computation easier.
+
+ 0 FFT_ENC-1
+ |-------------------------|
+
+ ----\ /----
+ \ /
+ \ / <- shifted version of window w[n]
+ \ /
+ \ /
+ -------
+
+ |---------| |---------|
+ NW/2 NW/2
+ */
+
+ for(i=0; i<FFT_ENC; i++) {
+ W[i].real = 0.0;
+ W[i].imag = 0.0;
+ }
+ for(i=0; i<NW/2; i++)
+ W[i].real = w[i+M/2];
+ for(i=FFT_ENC-NW/2,j=M/2-NW/2; i<FFT_ENC; i++,j++)
+ W[i].real = w[j];
+
+ four1(&W[-1].imag,FFT_ENC,-1); /* "Numerical Recipes in C" FFT */
+
+ /*
+ Re-arrange W[] to be symmetrical about FFT_ENC/2. Makes later
+ analysis convenient.
+
+ Before:
+
+
+ 0 FFT_ENC-1
+ |----------|---------|
+ __ _
+ \ /
+ \_______________/
+
+ After:
+
+ 0 FFT_ENC-1
+ |----------|---------|
+ ___
+ / \
+ ________/ \_______
+
+ */
+
+
+ for(i=0; i<FFT_ENC/2; i++) {
+ temp.real = W[i].real;
+ temp.imag = W[i].imag;
+ W[i].real = W[i+FFT_ENC/2].real;
+ W[i].imag = W[i+FFT_ENC/2].imag;
+ W[i+FFT_ENC/2].real = temp.real;
+ W[i+FFT_ENC/2].imag = temp.imag;
+ }
+
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: dft_speech
+ AUTHOR......: David Rowe
+ DATE CREATED: 27/5/94
+
+ Finds the DFT of the current speech input speech frame.
+
+\*---------------------------------------------------------------------------*/
+
+void dft_speech(COMP Sw[], float Sn[], float w[])
+{
+ int i;
+
+ for(i=0; i<FFT_ENC; i++) {
+ Sw[i].real = 0.0;
+ Sw[i].imag = 0.0;
+ }
+
+ /* Centre analysis window on time axis, we need to arrange input
+ to FFT this way to make FFT phases correct */
+
+ /* move 2nd half to start of FFT input vector */
+
+ for(i=0; i<NW/2; i++)
+ Sw[i].real = Sn[i+M/2]*w[i+M/2];
+
+ /* move 1st half to end of FFT input vector */
+
+ for(i=0; i<NW/2; i++)
+ Sw[FFT_ENC-NW/2+i].real = Sn[i+M/2-NW/2]*w[i+M/2-NW/2];
+
+ four1(&Sw[-1].imag,FFT_ENC,-1);
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: two_stage_pitch_refinement
+ AUTHOR......: David Rowe
+ DATE CREATED: 27/5/94
+
+ Refines the current pitch estimate using the harmonic sum pitch
+ estimation technique.
+
+\*---------------------------------------------------------------------------*/
+
+void two_stage_pitch_refinement(MODEL *model, COMP Sw[])
+{
+ float pmin,pmax,pstep; /* pitch refinment minimum, maximum and step */
+
+ /* Coarse refinement */
+
+ pmax = TWO_PI/model->Wo + 5;
+ pmin = TWO_PI/model->Wo - 5;
+ pstep = 1.0;
+ hs_pitch_refinement(model,Sw,pmin,pmax,pstep);
+
+ /* Fine refinement */
+
+ pmax = TWO_PI/model->Wo + 1;
+ pmin = TWO_PI/model->Wo - 1;
+ pstep = 0.25;
+ hs_pitch_refinement(model,Sw,pmin,pmax,pstep);
+
+ /* Limit range */
+
+ if (model->Wo < TWO_PI/P_MAX)
+ model->Wo = TWO_PI/P_MAX;
+ if (model->Wo > TWO_PI/P_MIN)
+ model->Wo = TWO_PI/P_MIN;
+
+ model->L = floor(PI/model->Wo);
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: hs_pitch_refinement
+ AUTHOR......: David Rowe
+ DATE CREATED: 27/5/94
+
+ Harmonic sum pitch refinement function.
+
+ pmin pitch search range minimum
+ pmax pitch search range maximum
+ step pitch search step size
+ model current pitch estimate in model.Wo
+
+ model refined pitch estimate in model.Wo
+
+\*---------------------------------------------------------------------------*/
+
+void hs_pitch_refinement(MODEL *model, COMP Sw[], float pmin, float pmax, float pstep)
+{
+ int m; /* loop variable */
+ int b; /* bin for current harmonic centre */
+ float E; /* energy for current pitch*/
+ float Wo; /* current "test" fundamental freq. */
+ float Wom; /* Wo that maximises E */
+ float Em; /* mamimum energy */
+ float r; /* number of rads/bin */
+ float p; /* current pitch */
+
+ /* Initialisation */
+
+ model->L = PI/model->Wo; /* use initial pitch est. for L */
+ Em = 0.0;
+ r = TWO_PI/FFT_ENC;
+
+ /* Determine harmonic sum for a range of Wo values */
+
+ for(p=pmin; p<=pmax; p+=pstep) {
+ E = 0.0;
+ Wo = TWO_PI/p;
+
+ /* Sum harmonic magnitudes */
+
+ for(m=1; m<=model->L; m++) {
+ b = floor(m*Wo/r + 0.5);
+ E += Sw[b].real*Sw[b].real + Sw[b].imag*Sw[b].imag;
+ }
+
+ /* Compare to see if this is a maximum */
+
+ if (E > Em) {
+ Em = E;
+ Wom = Wo;
+ }
+ }
+
+ model->Wo = Wom;
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: estimate_amplitudes
+ AUTHOR......: David Rowe
+ DATE CREATED: 27/5/94
+
+ Estimates the complex amplitudes of the harmonics.
+
+\*---------------------------------------------------------------------------*/
+
+void estimate_amplitudes(MODEL *model, COMP Sw[], COMP W[])
+{
+ int i,m; /* loop variables */
+ int am,bm; /* bounds of current harmonic */
+ int b; /* DFT bin of centre of current harmonic */
+ float den; /* denominator of amplitude expression */
+ float r; /* number of rads/bin */
+ int offset;
+ COMP Am;
+
+ r = TWO_PI/FFT_ENC;
+
+ for(m=1; m<=model->L; m++) {
+ den = 0.0;
+ am = floor((m - 0.5)*model->Wo/r + 0.5);
+ bm = floor((m + 0.5)*model->Wo/r + 0.5);
+ b = floor(m*model->Wo/r + 0.5);
+
+ /* Estimate ampltude of harmonic */
+
+ den = 0.0;
+ Am.real = Am.imag = 0.0;
+ for(i=am; i<bm; i++) {
+ den += Sw[i].real*Sw[i].real + Sw[i].imag*Sw[i].imag;
+ offset = i + FFT_ENC/2 - floor(m*model->Wo/r + 0.5);
+ Am.real += Sw[i].real*W[offset].real;
+ Am.imag += Sw[i].imag*W[offset].real;
+ }
+
+ model->A[m] = sqrt(den);
+
+ /* Estimate phase of harmonic */
+
+ model->phi[m] = atan2(Sw[b].imag,Sw[b].real);
+ }
+}
+
+/*---------------------------------------------------------------------------*\
+
+ est_voicing_mbe()
+
+ Returns the error of the MBE cost function for a fiven F0.
+
+ Note: I think a lot of the operations below can be simplified as
+ W[].imag = 0 and has been normalised such that den always equals 1.
+
+\*---------------------------------------------------------------------------*/
+
+float est_voicing_mbe(
+ MODEL *model,
+ COMP Sw[],
+ COMP W[],
+ float f0,
+ COMP Sw_[], /* DFT of all voiced synthesised signal for f0 */
+ /* useful for debugging/dump file */
+ int *voiced
+)
+{
+ int i,l,al,bl,m; /* loop variables */
+ COMP Am; /* amplitude sample for this band */
+ int offset; /* centers Hw[] about current harmonic */
+ float den; /* denominator of Am expression */
+ float error; /* accumulated error between originl and synthesised */
+ float Wo; /* current "test" fundamental freq. */
+ int L;
+ float sig, snr;
+
+ sig = 0.0;
+ for(l=1; l<=model->L/4; l++) {
+ sig += model->A[l]*model->A[l];
+ }
+
+ for(i=0; i<FFT_ENC; i++) {
+ Sw_[i].real = 0.0;
+ Sw_[i].imag = 0.0;
+ }
+
+ L = floor((FS/2.0)/f0);
+ Wo = f0*(TWO_PI/FS);
+
+ error = 0.0;
+
+ /* Just test across the harmonics in the first 1000 Hz (L/4) */
+
+ for(l=1; l<=L/4; l++) {
+ Am.real = 0.0;
+ Am.imag = 0.0;
+ den = 0.0;
+ al = ceil((l - 0.5)*Wo*FFT_ENC/TWO_PI);
+ bl = ceil((l + 0.5)*Wo*FFT_ENC/TWO_PI);
+
+ /* Estimate amplitude of harmonic assuming harmonic is totally voiced */
+
+ for(m=al; m<bl; m++) {
+ offset = FFT_ENC/2 + m - l*Wo*FFT_ENC/TWO_PI + 0.5;
+ Am.real += Sw[m].real*W[offset].real + Sw[m].imag*W[offset].imag;
+ Am.imag += Sw[m].imag*W[offset].real - Sw[m].real*W[offset].imag;
+ den += W[offset].real*W[offset].real + W[offset].imag*W[offset].imag;
+ }
+
+ Am.real = Am.real/den;
+ Am.imag = Am.imag/den;
+
+ /* Determine error between estimated harmonic and original */
+
+ for(m=al; m<bl; m++) {
+ offset = FFT_ENC/2 + m - l*Wo*FFT_ENC/TWO_PI + 0.5;
+ Sw_[m].real = Am.real*W[offset].real - Am.imag*W[offset].imag;
+ Sw_[m].imag = Am.real*W[offset].imag + Am.imag*W[offset].real;
+ error += (Sw[m].real - Sw_[m].real)*(Sw[m].real - Sw_[m].real);
+ error += (Sw[m].imag - Sw_[m].imag)*(Sw[m].imag - Sw_[m].imag);
+ }
+ }
+
+ snr = 10.0*log10(sig/error);
+ if (snr > V_THRESH)
+ *voiced = 1;
+ else
+ *voiced = 0;
+
+ return snr;
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: make_synthesis_window
+ AUTHOR......: David Rowe
+ DATE CREATED: 11/5/94
+
+ Init function that generates the trapezoidal (Parzen) sythesis window.
+
+\*---------------------------------------------------------------------------*/
+
+void make_synthesis_window(float Pn[])
+{
+ int i;
+ float win;
+
+ /* Generate Parzen window in time domain */
+
+ win = 0.0;
+ for(i=0; i<N/2-TW; i++)
+ Pn[i] = 0.0;
+ win = 0.0;
+ for(i=N/2-TW; i<N/2+TW; win+=1.0/(2*TW), i++ )
+ Pn[i] = win;
+ for(i=N/2+TW; i<3*N/2-TW; i++)
+ Pn[i] = 1.0;
+ win = 1.0;
+ for(i=3*N/2-TW; i<3*N/2+TW; win-=1.0/(2*TW), i++)
+ Pn[i] = win;
+ for(i=3*N/2+TW; i<2*N; i++)
+ Pn[i] = 0.0;
+}
+
+/*---------------------------------------------------------------------------*\
+
+ FUNCTION....: synthesise
+ AUTHOR......: David Rowe
+ DATE CREATED: 20/2/95
+
+ Synthesise a speech signal in the frequency domain from the
+ sinusodal model parameters. Uses overlap-add a triangular window to
+ smoothly interpolate betwen frames.
+
+\*---------------------------------------------------------------------------*/
+
+void synthesise(
+ float Sn_[], /* time domain synthesised signal */
+ MODEL *model, /* ptr to model parameters for this frame */
+ float Pn[], /* time domain Parzen window */
+ int shift /* used to handle transition frames */
+)
+{
+ int i,l,j,b; /* loop variables */
+ COMP Sw_[FFT_DEC]; /* DFT of synthesised signal */
+
+ if (shift) {
+ /* Update memories */
+
+ for(i=0; i<N-1; i++) {
+ Sn_[i] = Sn_[i+N];
+ }
+ Sn_[N-1] = 0.0;
+ }
+
+ for(i=0; i<FFT_DEC; i++) {
+ Sw_[i].real = 0.0;
+ Sw_[i].imag = 0.0;
+ }
+
+ /* Now set up frequency domain synthesised speech */
+
+ for(l=1; l<=model->L; l++) {
+ b = floor(l*model->Wo*FFT_DEC/TWO_PI + 0.5);
+ Sw_[b].real = model->A[l]*cos(model->phi[l]);
+ Sw_[b].imag = model->A[l]*sin(model->phi[l]);
+ Sw_[FFT_DEC-b].real = Sw_[b].real;
+ Sw_[FFT_DEC-b].imag = -Sw_[b].imag;
+ }
+
+ /* Perform inverse DFT */
+
+ four1(&Sw_[-1].imag,FFT_DEC,1);
+
+ /* Overlap add to previous samples */
+
+ for(i=0; i<N-1; i++) {
+ Sn_[i] += Sw_[FFT_DEC-N+1+i].real*Pn[i];
+ }
+
+ if (shift)
+ for(i=N-1,j=0; i<2*N; i++,j++)
+ Sn_[i] = Sw_[j].real*Pn[i];
+ else
+ for(i=N-1,j=0; i<2*N; i++,j++)
+ Sn_[i] += Sw_[j].real*Pn[i];
+}
+
AUTHOR......: David Rowe
DATE CREATED: 1/11/94
- Header file for Sinusoidal coder.
+ Header file for sinusoidal analysis and synthesis functions.
\*---------------------------------------------------------------------------*/
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
-/*---------------------------------------------------------------------------*\
-
- INCLUDES
-
-\*---------------------------------------------------------------------------*/
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-
-#include "defines.h" /* defines for sinusoidal coder */
-#include "globals.h" /* external globals */
-
-/*---------------------------------------------------------------------------*\
-
- FUNCTIONS
-
-\*---------------------------------------------------------------------------*/
-
-/* functions in refine.c */
-
-void dft_speech();
-void two_stage_pitch_refinement();
-void hs_pitch_refinement(float pmin, float pmax, float pstep);
-
-/* functions in spec.c */
-
-void estimate_amplitudes();
-void estimate_voicing();
-void estimate_voicing_av();
-float voicing(int lower, int upper);
-
-/* functions in four1.c */
-
-void four1();
-
-/* functions in synth.c */
-
-void synthesise_mixed();
-
-/* functions in initenc.c and initdec.c */
-
-void init_encoder(void);
-float make_window();
-void init_decoder();
-
-/* functions in gasdev.c */
+#ifndef __SINE__
+#define __SINE__
-float gasdev();
+void make_analysis_window(float w[], COMP W[]);
+void dft_speech(COMP Sw[], float Sn[], float w[]);
+void two_stage_pitch_refinement(MODEL *model, COMP Sw[]);
+void estimate_amplitudes(MODEL *model, COMP Sw[], COMP W[]);
+float est_voicing_mbe(MODEL *model, COMP Sw[], COMP W[], float f0, COMP Sw_[], int *voiced);
+void make_synthesis_window(float Pn[]);
+void synthesise(float Sn_[], MODEL *model, float Pn[], int shift);
+#endif
*/
#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
#include <string.h>
+#include <math.h>
+
+#include "defines.h"
#include "sine.h"
-#include "quantise.h"
#include "dump.h"
-#include "phase.h"
#include "lpc.h"
-#include "synth.h"
+#include "quantise.h"
+#include "phase.h"
#include "postfilter.h"
#include "interp.h"
-#include "nlp.h"
/*---------------------------------------------------------------------------*\
int main(int argc, char *argv[])
{
FILE *fmodel; /* file of model parameters from encoder */
- FILE *fout; /* output speech file */
- FILE *fin; /* input speech file */
- short buf[N]; /* input/output buffer */
- int i; /* loop variable */
- int length; /* number of frames so far */
+ FILE *fout; /* output speech file */
+ FILE *fin; /* input speech file */
+ short buf[N]; /* input/output buffer */
+ float Sn[M]; /* float input speech samples */
+ COMP Sw[FFT_ENC]; /* DFT of w[] */
+ float w[M]; /* time domain hamming window */
+ COMP W[FFT_ENC]; /* DFT of w[] */
+ MODEL model;
+ float Pn[2*N]; /* trapezoidal synthesis window */
+ float Sn_[2*N]; /* synthesised speech */
+ int i; /* loop variable */
+ int frames;
+ int length; /* number of frames so far */
char out_file[MAX_STR];
int arg;
int lpc_model, order;
int lsp, lsp_quantiser;
- float ak[LPC_MAX_ORDER+1];
+ float ak[LPC_MAX+1];
int dump;
- int phase, phase_model;
+ int phase0;
float ex_phase[MAX_AMP+1];
int voiced, voiced_1, voiced_2;
MODEL model_1, model_2, model_3, model_synth, model_a, model_b;
int transition, decimate;
+ for(i=0; i<2*N; i++)
+ Sn_[i] = 0;
+
voiced_1 = voiced_2 = 0;
model_1.Wo = TWO_PI/P_MIN;
model_1.L = floor(PI/model_1.Wo);
lsp = switch_present("--lsp",argc,argv);
lsp_quantiser = 0;
- /* phase_model 0: zero phase
- phase_model 1: 1st order polynomial */
- phase = switch_present("--phase",argc,argv);
- if (phase) {
- phase_model = atoi(argv[phase+1]);
- assert((phase_model == 0) || (phase_model == 1));
+ phase0 = switch_present("--phase0",argc,argv);
+ if (phase0) {
ex_phase[0] = 0;
}
/* Initialise ------------------------------------------------------------*/
- init_decoder();
- init_encoder();
- Nw = 220;
- make_window(Nw);
+ make_analysis_window(w,W);
+ make_synthesis_window(Pn);
quantise_init();
/* Main loop ------------------------------------------------------------*/
for(i=0; i<N; i++)
Sn[i+M-N] = buf[i];
dump_Sn(Sn);
- dft_speech(); dump_Sw(Sw);
+ dft_speech(Sw, Sn, w); dump_Sw(Sw);
dump_model(&model);
- /* optional phase modelling - make sure this happens before LPC
- modelling of {Am} as first order model fit doesn't work well
- with LPC Modelled {Am} (not sure why - investigate later) */
+ /* optional zero-phase modelling */
- if (phase) {
+ if (phase0) {
float Wn[M]; /* windowed speech samples */
- float Rk[PHASE_LPC_ORD+1]; /* autocorrelation coeffs */
- float ak_phase[PHASE_LPC_ORD+1];/* LPCs */
- COMP H[MAX_AMP]; /* LPC freq domain samples */
- float n_min;
- COMP min_Am;
+ float ak_phase[PHASE_LPC+1]; /* autocorrelation coeffs */
+ float Rk[PHASE_LPC+1]; /* autocorrelation coeffs */
+ COMP Sw_[FFT_ENC];
- dump_phase(&model.phi[0]);
+ dump_phase(&model.phi[0], model.L);
/* Determine LPCs for phase modelling. Note that we may also
find the LPCs as part of the {Am} modelling, this can
for(i=0; i<M; i++)
Wn[i] = Sn[i]*w[i];
- autocorrelate(Wn,Rk,M,PHASE_LPC_ORD);
- levinson_durbin(Rk,ak_phase,PHASE_LPC_ORD);
+ autocorrelate(Wn,Rk,M,PHASE_LPC);
+ levinson_durbin(Rk,ak_phase,PHASE_LPC);
if (lpc_model)
- assert(order == PHASE_LPC_ORD);
-
- dump_ak(ak_phase, PHASE_LPC_ORD);
-
- {
- COMP Sw_[FFT_ENC];
- float sig = 0.0;
- float noise;
- int m;
- float candidate_f0;
- float f0,best_f0;
- float e,e_min;
- int i;
- float f0_min, f0_max;
- float f0_start, f0_end;
-
- for(m=1; m<=model.L/4; m++) {
- sig += model.A[m]*model.A[m];
- }
- f0_min = (float)8000.0/P_MAX;
- f0_max = (float)8000.0/P_MIN;
- candidate_f0 = (4000.0/PI)*model.Wo;
- f0_start = candidate_f0-5;
- f0_end = candidate_f0+5;
- if (f0_start < f0_min) f0_start = f0_min;
- if (f0_end > f0_max) f0_end = f0_max;
- e_min = 1E32;
-
- for(f0=f0_start; f0<=f0_end; f0+= 1) {
- e = test_candidate_mbe(Sw, f0, Sw_);
- if (e < e_min) {
- e_min = e;
- best_f0 = f0;
- }
- }
- noise = test_candidate_mbe(Sw,best_f0, Sw_);
- snr = sig/noise;
- dump_Sw_(Sw_);
- }
-
- phase_model_first_order(ak_phase, H, &n_min, &min_Am, &voiced);
- if (snr > 4.0)
- voiced = 1;
- else
- voiced = 0;
+ assert(order == PHASE_LPC);
+
+ dump_ak(ak_phase, PHASE_LPC);
+ /* determine voicing */
+ snr = est_voicing_mbe(&model, Sw, W, (FS/TWO_PI)*model.Wo, Sw_, &voiced);
+ dump_Sw_(Sw_);
dump_snr(snr);
- if (phase_model == 0) {
- /* just to make sure we are not cheating - kill all phases */
- for(i=0; i<MAX_AMP; i++)
- model.phi[i] = 0;
- if (hand_snr) {
- fscanf(fsnr,"%f\n",&snr);
- voiced = snr > 2.0;
- }
- phase_synth_zero_order(voiced, H, ex_phase, voiced);
- }
- if (phase_model == 1) {
- phase_synth_first_order(voiced, H, n_min, min_Am, voiced);
- }
+ /* just to make sure we are not cheating - kill all phases */
- if (postfilt)
+ for(i=0; i<MAX_AMP; i++)
+ model.phi[i] = 0;
+ if (hand_snr) {
+ fscanf(fsnr,"%f\n",&snr);
+ voiced = snr > 2.0;
+ }
+ phase_synth_zero_order(&model, ak_phase, voiced, ex_phase);
+
+ if (postfilt)
postfilter(&model, voiced, &bg_est);
-
- //dump_phase_(&model.phi[0]);
}
/* optional LPC model amplitudes */
if (lpc_model) {
- snr = lpc_model_amplitudes(Sn, &model, order, lsp, ak);
+ snr = lpc_model_amplitudes(Sn, w, &model, order, lsp, ak);
sum_snr += snr;
dump_quantised_model(&model);
}
if (fout != NULL) {
if (transition) {
- synthesise_mixed(Pn,&model_a,Sn_,1);
- synthesise_mixed(Pn,&model_b,Sn_,0);
+ synthesise(Sn_,&model_a,Pn,1);
+ synthesise(Sn_,&model_b,Pn,0);
}
else {
- synthesise_mixed(Pn,&model,Sn_,1);
+ synthesise(Sn_,&model,Pn,1);
}
/* Save output speech to disk */
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
-#include "sine.h"
-#include "dump.h"
+#include <stdio.h>
+#include <stdlib.h>
#include <string.h>
+#include "defines.h"
+#include "dump.h"
+#include "sine.h"
/*---------------------------------------------------------------------------*\
int main(int argc, char *argv[])
{
- FILE *fin; /* input speech sample file */
- FILE *fmodel; /* output file of model parameters */
- FILE *fp; /* input text file containing pitch estimates */
- short buf[N]; /* input speech sample buffer */
- int length; /* number of frames to process */
- float pitch; /* current pitch estimate from external pitch file */
- int i; /* loop variable */
- FILE *fref; /* optional output file with refined pitch estimate */
- int arg;
- int dump;
- int frames;
+ FILE *fin; /* input speech sample file */
+ FILE *fmodel; /* output file of model parameters */
+ FILE *fp; /* input text file containing pitch estimates */
+ short buf[N]; /* input speech sample buffer */
+ float Sn[M]; /* float input speech samples */
+ COMP Sw[FFT_ENC]; /* DFT of Sn[] */
+ float w[M]; /* time domain hamming window */
+ COMP W[FFT_ENC]; /* DFT of w[] */
+ MODEL model;
+ int length; /* number of frames to process */
+ float pitch; /* current pitch estimate from external pitch file */
+ int i; /* loop variable */
+ FILE *fref; /* optional output file with refined pitch estimate */
+ int arg;
+ int dump;
+ int frames;
if (argc < 5) {
printf("usage: sinenc InputFile ModelFile Frames PitchFile\n");
if (dump)
dump_on(argv[dump+1]);
- if ((arg == switch_present("--ref",argc,argv))) {
+ if ((arg = switch_present("--ref",argc,argv))) {
if ((fref = fopen(argv[arg+1],"wt")) == NULL) {
printf("Error opening output pitch refinement file: %s\n",argv[5]);
exit(1);
else
fref = NULL;
- init_encoder();
- make_window(NW);
+ /* Initialise sample buffer memories to stop divide by zero errors
+ and zero energy frames at the start of simulation */
+
+ for(i=0; i<M; i++)
+ Sn[i] = 1.0;
+ make_analysis_window(w, W);
+
/* Main loop ------------------------------------------------------------*/
frames = 0;
/* estimate and model parameters */
- dft_speech();
- two_stage_pitch_refinement();
- estimate_amplitudes();
- dump_Sn(Sn); dump_Sw(Sw); dump_Sw_(Sw_); dump_model(&model);
+ dft_speech(Sw, Sn, w);
+ two_stage_pitch_refinement(&model, Sw);
+ estimate_amplitudes(&model, Sw, W);
+ dump_Sn(Sn); dump_Sw(Sw); dump_model(&model);
/* save model parameters */
AUTHOR......: David Rowe
DATE CREATED: 27/5/94
- Functions for estimating the complex amplitude of harmonics.
+ Functions for estimating the sinusoidal model parameters.
\*---------------------------------------------------------------------------*/
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
-#include "sine.h"
+#include "defines.h"
+#include "spec.h"
/*---------------------------------------------------------------------------*\
AUTHOR......: David Rowe
DATE CREATED: 27/5/94
- Estimates the complex amplitudes of the harmonics. Also generates
- all voiced synthetic spectrum for later voicing estimation.
-
- INPUT.......: global float Sw[] DFT of speech
- global MODEL model contains parameters L and Wo
-
- OUTPUT......: global float Sw_[] DFT of all voiced synthesised speech
- global MODEL model contains parameters A[] and phi[]
+ Estimates the complex amplitudes of the harmonics.
\*---------------------------------------------------------------------------*/
-void estimate_amplitudes()
+void estimate_amplitudes(MODEL *model, float Sw[])
{
- int i,m; /* loop variables */
- int am,bm; /* bounds of current harmonic */
- int b; /* DFT bin of centre of current harmonic */
+ int i,m; /* loop variables */
+ int am,bm; /* bounds of current harmonic */
+ int b; /* DFT bin of centre of current harmonic */
float den; /* denominator of amplitude expression */
float r; /* number of rads/bin */
float E;
/* Estimate phase of harmonic */
model.phi[m] = atan2(Sw[b].imag,Sw[b].real);
-
- #ifdef MBE_VOICING_NEEDED
- /* construct all voiced model spectrum and estimate voicing using MBE model */
-
- E = 0.0;
- for(i=am; i<bm; i++) {
- offset = FFT_ENC/2 + i - floor(m*model.Wo/r + 0.5);
- Sw_[i].real = Am.real*W[offset].real;
- Sw_[i].imag = Am.imag*W[offset].real;
- E = pow(Sw[i].real - Sw_[i].real, 2.0) + pow(Sw[i].imag - Sw_[i].imag, 2.0);
- }
- model.v[m] = E/den;
- #endif
- }
}
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
-#include "sine.h"
+#include "defines.h"
+#include "synth.h"
-void synthesise_mixed(
- float Pn[], /* time domain Parzen window */
+void synthesise(
+ float Pn[], /* time domain Parzen window */
MODEL *model, /* ptr to model parameters for this frame */
- float Sn_[], /* time domain synthesised signal */
- int shift /* if non-zero update memories */
+ float Sn_[], /* time domain synthesised signal */
+ int shift /* used to handle transition frames */
)
{
- int i,l,j,b; /* loop variables */
- COMP Sw_[FFT_DEC]; /* DFT of synthesised signal */
+ int i,l,j,b; /* loop variables */
+ COMP Sw_[FFT_DEC]; /* DFT of synthesised signal */
- if (shift) {
- /* Update memories */
+ if (shift) {
+ /* Update memories */
- for(i=0; i<N-1; i++) {
- Sn_[i] = Sn_[i+N];
- }
- Sn_[N-1] = 0.0;
- }
+ for(i=0; i<N-1; i++) {
+ Sn_[i] = Sn_[i+N];
+ }
+ Sn_[N-1] = 0.0;
+ }
- for(i=0; i<FFT_DEC; i++) {
- Sw_[i].real = 0.0;
- Sw_[i].imag = 0.0;
- }
+ for(i=0; i<FFT_DEC; i++) {
+ Sw_[i].real = 0.0;
+ Sw_[i].imag = 0.0;
+ }
- /* Now set up frequency domain synthesised speech */
+ /* Now set up frequency domain synthesised speech */
- for(l=1; l<=model->L; l++) {
- b = floor(l*model->Wo*FFT_DEC/TWO_PI + 0.5);
- Sw_[b].real = model->A[l]*cos(model->phi[l]);
- Sw_[b].imag = model->A[l]*sin(model->phi[l]);
- Sw_[FFT_DEC-b].real = Sw_[b].real;
- Sw_[FFT_DEC-b].imag = -Sw_[b].imag;
- }
+ for(l=1; l<=model->L; l++) {
+ b = floor(l*model->Wo*FFT_DEC/TWO_PI + 0.5);
+ Sw_[b].real = model->A[l]*cos(model->phi[l]);
+ Sw_[b].imag = model->A[l]*sin(model->phi[l]);
+ Sw_[FFT_DEC-b].real = Sw_[b].real;
+ Sw_[FFT_DEC-b].imag = -Sw_[b].imag;
+ }
- /* Perform inverse DFT */
+ /* Perform inverse DFT */
- four1(&Sw_[-1].imag,FFT_DEC,1);
+ four1(&Sw_[-1].imag,FFT_DEC,1);
- /* Overlap add to previous samples */
+ /* Overlap add to previous samples */
- for(i=0; i<N-1; i++) {
- Sn_[i] += Sw_[FFT_DEC-N+1+i].real*Pn[i];
- }
+ for(i=0; i<N-1; i++) {
+ Sn_[i] += Sw_[FFT_DEC-N+1+i].real*Pn[i];
+ }
- if (shift)
- for(i=N-1,j=0; i<2*N; i++,j++)
- Sn_[i] = Sw_[j].real*Pn[i];
- else
- for(i=N-1,j=0; i<2*N; i++,j++)
- Sn_[i] += Sw_[j].real*Pn[i];
-
-}
-
-/*---------------------------------------------------------------------------*\
-
- synthesise_continuous_phase()
-
- This version adjusts the frequency of each harmonic slightly to
- ensure a continuous phase track from the previous frame. Used with
- the zero phase model, when original phases are not available.
-
- At sample n=0 of this frame, we assume the phase of harmonic m is
- set to phi_prev[m]. We want the final phase at sample N to be
- phi[m]. This means the phase track must start at phi_prev[m],
- rotate several times based on mWo, then end up at phase phi[m].
-
- To ensure that the phase track arrives at phi[m] by sample N we add
- a small frequency offset by slightly shifting the frequency of each
- harmonic.
-
- The continuous phase track model is only useful for voiced speech.
- In fact, for unvoiced speech we desire a rough, discontinuous phase
- track. So in unvoiced frames or in cases where the fundamental
- frequency varies by more that 20%, we don't add the small frequency
- offset.
-
- Result: when tested was no difference in output speech quality. The
- partial unvoiced sound when using zero phase model was found to be
- due mis-alignment of the LPC analysis window and accidental addition
- of a random phase component. So we are sticking with synthesise_mixed()
- above for now. I am leaving this function here for the moment as it
- might be useful one day.
-
-\*---------------------------------------------------------------------------*/
-
-void synthesise_continuous_phase(
- float Pn[], /* time domain Parzen window */
- MODEL *model, /* ptr to model parameters for this frame */
- float Sn_[], /* time domain synthesised signal */
- int voiced, /* non-zero if frame is voiced */
- float *Wo_prev, /* previous frames Wo */
- float phi_prev[] /* previous frames phases */
-)
-{
- int i,l,j; /* loop variables */
- COMP Sw_[FFT_DEC]; /* DFT of synthesised signal */
- int b[MAX_AMP]; /* DFT bin of each harmonic */
- float delta_w; /* frequency offset required */
-
- /* Update memories */
-
- for(i=0; i<N-1; i++) {
- Sn_[i] = Sn_[i+N];
- }
- Sn_[N-1] = 0.0;
-
- for(i=0; i<FFT_DEC; i++) {
- Sw_[i].real = 0.0;
- Sw_[i].imag = 0.0;
- }
-
- if (!voiced || (fabs(*Wo_prev - model->Wo) > 0.2*model->Wo)) {
- //printf("disc voiced = %d\n", voiced);
- //printf("%f %f\n", fabs(*Wo_prev - model->Wo), 0.2*model->Wo);
- /* discontinous phase tracks: no phase adjustment of frequency
- as we want discontinuous phase tracks */
-
- for(l=1; l<=model->L; l++)
- b[l] = floor(l*model->Wo*FFT_DEC/TWO_PI + 0.5);
- }
- else {
- /* continous phase tracks: determine frequency of each harmonic
- to ensure smooth phase track at the centre of next synthesis
- frame */
-
- for(l=1; l<=model->L; l++) {
- //printf("Wo_prev = %f Wo = %f\n", *Wo_prev, model->Wo);
- delta_w = (model->phi[l] - l*N*(*Wo_prev + model->Wo)/2.0 - phi_prev[l]);
- delta_w -= TWO_PI*floor(delta_w/TWO_PI + 0.5);
- delta_w /= N;
- b[l] = floor((l*model->Wo+delta_w)*FFT_DEC/TWO_PI + 0.5);
- //printf("delta_w = %f b[%d] = %d\n", delta_w,l,b[l]);
- }
- }
-
- /* update memories for phase tracking */
-
- *Wo_prev = model->Wo;
- for(l=1; l<=model->L; l++)
- phi_prev[l] = model->phi[l];
-
- /* Now set up frequency domain synthesised speech */
-
- for(l=1; l<=model->L; l++) {
- Sw_[b[l]].real = model->A[l]*cos(model->phi[l]);
- Sw_[b[l]].imag = model->A[l]*sin(model->phi[l]);
- Sw_[FFT_DEC-b[l]].real = Sw_[b[l]].real;
- Sw_[FFT_DEC-b[l]].imag = -Sw_[b[l]].imag;
- }
-
- /* Perform inverse DFT */
-
- four1(&Sw_[-1].imag,FFT_DEC,1);
-
- /* Overlap add to previous samples */
-
- for(i=0; i<N-1; i++) {
- Sn_[i] += Sw_[FFT_DEC-N+1+i].real*Pn[i];
- }
- for(i=N-1,j=0; i<2*N; i++,j++)
- Sn_[i] = Sw_[j].real*Pn[i];
+ if (shift)
+ for(i=N-1,j=0; i<2*N; i++,j++)
+ Sn_[i] = Sw_[j].real*Pn[i];
+ else
+ for(i=N-1,j=0; i<2*N; i++,j++)
+ Sn_[i] += Sw_[j].real*Pn[i];
}
#include "sine.h"
-void synthesise_mixed(float Pn[], MODEL *model, float Sn_[], int shift);
-void synthesise_continuous_phase(float Pn[], MODEL *model, float Sn_[],
- int voiced, float *Wo_prev, float phi_prev[]);
+void synthesise(float Pn[], MODEL *model, float Sn_[], int shift);
+
#endif
--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ FILE........: window.c
+ AUTHOR......: David Rowe
+ DATE CREATED: 11/5/94
+
+ Generates the time domain analysis window and it's DFT.
+
+\*---------------------------------------------------------------------------*/
+
+/*
+ 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 <math.h>
+#include "defines.h"
+#include "window.h"
+
+float make_window(
+ float w[], /* time domain analysis window */
+ COMP W[] /* w[] in frequency domain */
+)
+{
+ float m;
+ COMP temp;
+ int i,j;
+
+ /*
+ Generate Hamming window centered on M-sample pitch analysis window
+
+ 0 M/2 M-1
+ |-------------|-------------|
+ |-------|-------|
+ NW samples
+
+ All our analysis/synthsis is centred on the M/2 sample.
+ */
+
+ m = 0.0;
+ for(i=0; i<M/2-NW/2; i++)
+ w[i] = 0.0;
+ for(i=M/2-NW/2,j=0; i<M/2+NW/2; i++,j++) {
+ w[i] = 0.5 - 0.5*cos(TWO_PI*j/(NW-1));
+ m += w[i]*w[i];
+ }
+ for(i=M/2+NW/2; i<M; i++)
+ w[i] = 0.0;
+
+ /* Normalise - makes freq domain amplitude estimation straight
+ forward */
+
+ m = 1.0/sqrt(m*FFT_ENC);
+ for(i=0; i<M; i++) {
+ w[i] *= m;
+ }
+
+ /*
+ Generate DFT of analysis window, used for later processing. Note
+ we modulo FFT_ENC shift the time domain window w[], this makes the
+ imaginary part of the DFT W[] equal to zero as the shifted w[] is
+ even about the n=0 time axis if NW is odd. Having the imag part
+ of the DFT W[] makes computation easier.
+
+ 0 FFT_ENC-1
+ |-------------------------|
+
+ ----\ /----
+ \ /
+ \ / <- shifted version of window w[n]
+ \ /
+ \ /
+ -------
+
+ |---------| |---------|
+ NW/2 NW/2
+ */
+
+ for(i=0; i<FFT_ENC; i++) {
+ W[i].real = 0.0;
+ W[i].imag = 0.0;
+ }
+ for(i=0; i<NW/2; i++)
+ W[i].real = w[i+M/2];
+ for(i=FFT_ENC-NW/2,j=M/2-NW/2; i<FFT_ENC; i++,j++)
+ W[i].real = w[j];
+
+ four1(&W[-1].imag,FFT_ENC,-1); /* "Numerical Recipes in C" FFT */
+
+ /*
+ Re-arrange W[] to be symmetrical about FFT_ENC/2. Makes later
+ analysis convenient.
+
+ Before:
+
+
+ 0 FFT_ENC-1
+ |----------|---------|
+ __ _
+ \ /
+ \_______________/
+
+ After:
+
+ 0 FFT_ENC-1
+ |----------|---------|
+ ___
+ / \
+ ________/ \_______
+
+ */
+
+
+ for(i=0; i<FFT_ENC/2; i++) {
+ temp.real = W[i].real;
+ temp.imag = W[i].imag;
+ W[i].real = W[i+FFT_ENC/2].real;
+ W[i].imag = W[i+FFT_ENC/2].imag;
+ W[i+FFT_ENC/2].real = temp.real;
+ W[i+FFT_ENC/2].imag = temp.imag;
+ }
+
+ return(m);
+}
+
--- /dev/null
+/*---------------------------------------------------------------------------*\
+
+ FILE........: window.h
+ AUTHOR......: David Rowe
+ DATE CREATED: 16/8/2010
+
+ Generates the time domain analysis window and it's DFT.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef __WINDOW__
+#define __WINDOW__
+
+float make_window(float w[], COMP W[]);
+
+#endif
projects / freetel-svn-tracking.git / commitdiff