for( i = 0, re = f, im = f + fftsz; i <= m2; i++ ){            for( j = 0; j < fftsz; j++ )	        *(re++) =  cos( ww[j] * i ) * dw[j];            for( j = 0; j < fftsz; j++ )	        *(im++) = -sin( ww[j] * i ) * dw[j];	    re -= fftsz;            im -= fftsz;            ifft( re, im, fftsz );            for( j = 1; j <= m1; j++ )	      re[j] += re[fftsz-j];            re += fftsz2;            im += fftsz2;        }        free( ww );        free( dw );        /*-------  copy "f" to "g" ----------*/                    for( i = 0, next = f, pf = f, pg = g; i <= m2; i++ ){            for( j = 0; j <= m1; j++ )                *(pg++) = *(pf++);            next += fftsz2;            pf = next;        }        free( f );        flag_g = 1;        for( j = 1; j <= m1; j++ )            g[j] *= 0.5;        for( i = 1; i <= m2; i++ )            g[i*(m1+1)] *= 2.0;            }    for( i = 0, pg = g; i <= m2; i++ )        for( j = 0, c2[i] = 0.; j <= m1; j++ )            c2[i] += *(pg++) * c1[j];}/***************************************************************  No.3  ifreqt2    static : *h, size3   Inverse Frequency Transformation	void	ifreqt2(c1, m1, c2, m2, fftsz, a, t)	double	*c1   : minimum phase sequence	int	m1    : order of minimum phase sequence	double	*c2   : warped sequence	int	m2    : order of warped sequence        int	fftsz : ifft size	double	a     : all-pass constant        double	t     : emphasized frequency t * pi(rad)***************************************************************/void ifreqt2(c1, m1, c2, m2, fftsz, a, t)double *c1, *c2, a, t;int m1, m2, fftsz;{    register int 	i, j;    double		w, b, *ww, *f,                        *re, *im, *pl, *pr, *plnxt, *prnxt,                        *pf, *ph, *next, warp(), derivw();    int                 size_h, size_f, fftsz2, m12, m11;    static double	*h = NULL;    static int		size3, flag_h = 1;    b = M_2PI / (double)fftsz;    size_h = ( m2 + 1 ) * ( m1 + 1 );    if( h == NULL ){        flag_h = 0;        size3 = size_h;        h = dgetmem( size3 );    }    else if( size_h != size3 ){        free(h);        flag_h = 0;                size3 = size_h;        h = dgetmem( size3 );    }    /*-------  if "h" is not defined  ----------*/        if( flag_h == 0 ){        ww = dgetmem( fftsz );        for( j = 0, w = 0.; j < fftsz; j++, w+=b )            ww[j] = warp( w, a, t );        fftsz2 = fftsz + fftsz;       /* size of (re + im) */	m12 = m1 + m1 + 1;        size_f = m12 * fftsz2;      /* size of array "f" */        f = dgetmem( size_f );                for( i = -m1, re = f, im = f + fftsz; i <= m1; i++ ){	  for( j = 0; j < fftsz; j++ )	    *(re++) =  cos( ww[j] * i );	  	  for( j = 0; j < fftsz; j++ )	    *(im++) = -sin( ww[j] * i );	  re -= fftsz;	  im -= fftsz;	  ifft( re, im, fftsz );	  re += fftsz2;	  im += fftsz2;        }        free( ww );        /*------- b'(n,m)=b(n,m)+b(n,-m) ----------*/            	pl = f;	pr = f + ( m12 - 1 ) * fftsz2;	for( i = 0, plnxt = pl, prnxt = pr; i < m1; i++ ){	  plnxt += fftsz2;	  prnxt -= fftsz2;	  for( j = 0; j <= m2; j++ )	    *(pr++) += *(pl++);	  pl = plnxt;	  pr = prnxt;	}        /*-------  copy "f" to "h" ----------*/            	m11 = m1 + 1;	pf = f + m1 * fftsz2;	for( j = 0, next = pf;  j <= m1; j++ ){	  next += fftsz2;	  for( i = 0; i <= m2; i++ )	    h[m11*i+j] = *(pf++);	  pf = next;	}        free( f );        flag_h = 1;        for( j = 1; j <= m1; j++ )            h[j] *= 0.5;        for( i = 1; i <= m2; i++ )            h[i*m11] *= 2.0;    }    for( i = 0, ph = h; i <= m2; i++ )      for( j = 0, c2[i] = 0.; j <= m1; j++ )	c2[i] += *(ph++) * c1[j];}/***************************************************************  No.4  frqtr2    static : *k, size4    Frequency Transformation for Calculating Coefficients	void	frqtr2(c1, m1, c2, m2, fftsz, a, t)	double	*c1   : minimum phase sequence	int	m1    : order of minimum phase sequence	double	*c2   : warped sequence	int	m2    : order of warped sequence	int	fftsz  : frame length (fft size)	double	a     : all-pass constant	double	t     : emphasized frequency***************************************************************/void frqtr2(c1, m1, c2, m2, fftsz, a, t)double *c1, *c2, a, t;int m1, m2, fftsz;{    register int 	i, j;    double		w, b, *ww, *f, *tc2,                        *re, *im, *pf, *pk, *next, warp();    int                 size_k, size_f, fftsz2;    static double	*k = NULL;    static int		size4, flag_k = 1;    b = M_2PI / (double)fftsz;    size_k = ( m2 + 1 ) * ( m1 + 1 );    if( k == NULL ){        flag_k = 0;        size4 = size_k;        k = dgetmem( size4 );    }    else if( size_k != size4 ){        free(k);        flag_k = 0;                size4 = size_k;        k = dgetmem( size4 );    }    /*-------  if "k" is not defined  ----------*/        if( flag_k == 0 ){        ww = dgetmem( fftsz );        for( j = 0, w = 0.; j < fftsz; j++, w+=b )            ww[j] = warp( w, a, t );        fftsz2 = fftsz + fftsz;         /* size of (re + im) */        size_f = ( m2 + 1 ) * fftsz2; /* size of array "f" */        f = dgetmem( size_f );                for( i = 0, re = f, im = f + fftsz; i <= m2; i++ ){            for( j = 0; j < fftsz; j++ )	      *(re++) =  cos( ww[j] * i );            for( j = 0; j < fftsz; j++ )	      *(im++) = -sin( ww[j] * i );	    re -= fftsz;            im -= fftsz;            ifft( re, im, fftsz );            for( j = 1; j <= m1; j++ )	      re[j] += re[fftsz-j];	                re += fftsz2;            im += fftsz2;        }        free( ww );        /*-------  copy "f" to "k" ----------*/                    for( i = 0, next = f, pf = f, pk = k; i <= m2; i++ ){            for( j = 0; j <= m1; j++ )                *(pk++) = *(pf++);            next += fftsz2;            pf = next;        }        free( f );        flag_k = 1;    }    tc2 = dgetmem( m2 + 1 );   /*  tmp of c2  */        for( i = 0, pk = k; i <= m2; i++ )        for( j = 0, tc2[i] = 0.; j <= m1; j++ )            tc2[i] += *(pk++) * c1[j];    movem(tc2, c2, sizeof(*c2), m2+1);    free(tc2);}/***************************************************************  Warping Function and Its Derivative    double   warp(w, a, t)  &  derivw(w, a, t)  double  w     : frequency  double  a     : all-pass constant  double  t     : emphasized frequency  ***************************************************************/double warp(w, a, t)double w, a, t;{    double              ww, x, y;    x = w-t;    y = w+t;        ww = w + atan2( (a * sin(x)), (1. - a * cos(x)) )           + atan2( (a * sin(y)), (1. - a * cos(y)) );    return(ww);}/*============================================================*/double derivw(w, a, t)double w, a, t;{    double dw, x, y, a2, aa;     x = w-t;    y = w+t;    a2 = a+a;    aa = a*a;    dw = 1. + ( a * cos(x) - aa )/( 1. - a2 * cos(x) + aa )            + ( a * cos(y) - aa )/( 1. - a2 * cos(y) + aa );    return(dw);}