/* Sound_and_Wavelet.c
 *
 * Copyright (C) 1992-2002 David Weenink & Paul Boersma
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * 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.
 */

/*
 * pb 1994/09/14
 * pb 2002/07/16 GPL
 */

#include "Sound.h"
#include "Wavelet.h"
#include "Sound_and_Wavelet.h"

static int NUMdwtset (int ncof, float **cc, float **cr);
static int NUMdwt (int ncof, float a [], long n, int isign);

Wavelet Sound_to_Wavelet (I, int ncof) {
	iam (Sound);
	Wavelet thee = NULL;
	long n = 1;
	while (n < my nx) n *= 2;
	thee = Wavelet_create (n, 0.5 / my dx, ncof);
	if (! thee) return NULL;
	thy highestFrequency = 0.5 / my dx;
	NUMfvector_copyElements (my z [1], thy a, 1, my nx);
	if (! NUMdwt (thy ncof, thy a, n, 1)) {
		forget (thee);
		return Melder_errorp ("(Sound_to_Wavelet:) Cannot create Wavelet transform.");
	}
	return thee;
} 
 
Sound Wavelet_to_Sound (I) {
	iam (Wavelet);
	Sound thee;
	thee = Sound_createSimple (my n / (2.0 * my highestFrequency), 2.0 * my highestFrequency);
	if (! thee) return NULL;
	NUMfvector_copyElements (my a, thy z [1], 1, my n);
	if (! NUMdwt (my ncof, thy z [1], my n, -1))
	{
		forget (thee);
		return Melder_errorp ("(Wavelet_to_Sound:) Cannot create Wavelet transform.");
	}
	return thee;
}

/* Initializing routine for NUMdwt, here implementing the Daubechies wavelet
 * filters with 4,6,8,10,12,14,16,18 and 20 coefficients, as selected by the 
 * input value ncof. The cc and cr point to the array of coefficients for the 
 * forward and reverse transform, respectively. 
 * The coefficients can be found in table 6.1 page 195 in: 
 *    Daubechies, I. (1992), "Ten lectures on wavelets",CBMS-NSF regional
 *    conference series in applied mathematics 61.
 */  
static int NUMdwtset(int ncof, float **cc, float **cr) {
	 static float cn[117] = { 
			 /* ncof == 4 */    0.0,
			0.4829629131445341, 0.8365163037378079, 0.2241438680420134,
		  -0.1294095225512604,  
			/* ncof ==  6 */    0.0,
			0.3326705529500825, 0.8068915093110924, 0.4598775021184914,
		  -0.1350110200102546,-0.0854412738820267, 0.0352262918857095,
			/* ncof ==  8 */    0.0,
			0.2303778133088964, 0.7148465705529154, 0.6308807679398587,
		  -0.0279837694168599,-0.1870348117190931, 0.0308413818355607,
			0.0328830116668852,-0.0105974017850690,
			/* ncof == 10 */    0.0,
			0.1601023979741929, 0.6038292697971895, 0.7243085284377726,
			0.1384281459013203,-0.2422948870663823,-0.0322448695846381,
			0.0775714938400459,-0.0062414902127983,-0.0125807519990820,
			0.0033357252854738,  
			/* ncof == 12 */    0.0,
			0.1115407433501095, 0.4946238903984533, 0.7511339080210959,
			0.3152503517091982,-0.2262646939654400,-0.1297668675672625,
			0.0975016055873225, 0.0275228655303053,-0.0315820393174862,
			0.0005538422011614, 0.0047772575109455,-0.0010773010853085,
			/* ncof == 14 */    0.0,
			0.0778520540850037, 0.3965393194818912, 0.7291320908461957,
			0.4697822874051889,-0.1439060039285212,-0.2240361849938412,
			0.0713092192668272, 0.0806126091510774,-0.0380299369350104,
		  -0.0165745416306655, 0.0125509985560986, 0.0004295779729214,
		  -0.0018016407040473, 0.0003537137999745,
			/* ncof == 16 */    0.0,
			0.0544158422431072, 0.3128715909143166, 0.6756307362973195,
			0.5853546836542159,-0.0158291052563823,-0.2840155429615824,
			0.0004724845739124, 0.1287474266204893,-0.0173693010018090,
		  -0.0440882539307971, 0.0139810279174001, 0.0087460940474065,
		  -0.0048703529934520,-0.0003917403733770, 0.0006754494064506,
		  -0.0001174767841248,
			/* ncof == 18 */    0.0,
			0.0380779473638778, 0.2438346746125858, 0.6048231236900955,
			0.6572880780512736, 0.1331973858249883,-0.2932737832791663,
		  -0.0968407832229492, 0.1485407493381256, 0.0307256814793385,
		  -0.0676328290613279, 0.0002509471148340, 0.0223616621236798,
		  -0.0047232047577518,-0.0042815036824635, 0.0018476468830563,
			0.0002303857635232,-0.0002519631889427, 0.0000393473203163,
			/* ncof == 20 */    0.0,
			0.0266700579005473, 0.1881768000776347, 0.5272011889315757,
			0.6884590394534363, 0.2811723436605715,-0.2498464243271598,
		  -0.1959462743772862, 0.1273693403357541, 0.0930573646035547,
		  -0.0713941471663501,-0.0294575368218399, 0.0332126740593612,
			0.0036065535669870,-0.0107331754833007, 0.0013953517470688,
			0.0019924052951925,-0.0006858566949564,-0.0001164668551285,
			0.0000935886703202,-0.0000132642028945 };
	 static float cnr[5+7+9+11+13+15+17+19+21];
	 int k, index;
	 float sig = -1.0;

	 if( ncof < 4 || ncof > 20 || ncof % 2 == 1) return 0;
	 index = ncof*ncof/4 - 4;
	 *cc = cn+index;
	 *cr = cnr+index;
	 for ( k=1; k <= ncof; k++ ) {
		  (*cr)[ncof+1-k] = sig * (*cc)[k];
		  sig = -sig;
	 }
	 return 1;
}

/* One-dimensional Discrete wavelet transform.
 * Function:
 *     This routine implements the pyramid algorithm, replacing a[1..na] by
 *     its wavelet transform (for isign=1), or performing the inverse operation
 *     (for isign=-1). Note that na must be an integer power of 2.
 *     ncof determines the number of filter coefficients used in the transform.
 */
static int NUMdwt(int ncof, float a[], long na, int isign) {
	 float ai, ai1, *wksp;
	 unsigned long i, ii, j, jf, jr, k, n1, ni, nj, nh, n, nmod;
	 int ioff = -(ncof >> 1); 
	 int joff = ioff;
	 float *cc, *cr;

	 /* get the Daubechies filter coefficients */
	 if ( na < 4 || ! NUMdwtset(ncof, &cc, &cr ) ) return 0;

	 wksp = NUMfvector(1,na);
	 if( !wksp ) return 0;
	 if (isign >= 0) {
	 /* filter: start at largest hierarchy, work towards smallest */
		for( n=na; n >= 4; n>>=1 ) {
			nmod = ncof*n;
			n1 = n-1;
			nh = n >> 1;
			for ( j=1; j <= n; j++ ) wksp[j] = 0.0;
			for ( ii=1, i=1; i <= n; i+=2, ii++ ) {
					 ni = i + nmod + ioff;
					 nj = i + nmod + joff;
					 for ( k=1; k <= ncof; k++ )
					 {
						  jf = n1 & (ni+k);
						  jr = n1 & (nj+k);
						  wksp[ii] +=  cc[k] * a[jf+1];
						  wksp[ii+nh] +=  cr[k] * a[jr+1];
					 }
				}
				for ( j=1; j <= n; j++ ) a[j] = wksp[j];
		  }
	 } else {
	 /* transpose filter: start at smallest hierarchy, work towards largest */
		  for( n=4; n <= na; n<<=1 ) {
			nmod = ncof*n;
			n1 = n-1;
			nh = n >> 1;
			for ( j=1; j <= n; j++ ) wksp[j] = 0.0;
			for ( ii=1, i=1; i <= n; i+=2, ii++ ) {
				ai = a[ii];
				ai1 = a[ii+nh];
				ni = i+nmod+ ioff;
				nj = i+nmod+ joff;
				for (k=1;k<= ncof;k++) {
					jf = (n1 & (ni+k))+1;
					jr = (n1 & (nj+k))+1;
					wksp[jf] +=  cc[k]*ai;
					wksp[jr] +=  cr[k]*ai1;
				}
			}
			for ( j=1; j <= n; j++ ) a[j] = wksp[j];
		  }
	 }
	 NUMfvector_free(wksp,1);

	 return 1;
}

/* End of file Sound_and_Wavelet.c */