bsynz.c,v

这是LPC-10压缩算法的源代码,愿共享之.也希望能赐与MELP方面的算法源码.

    /* Function Body */
    ipo = &(st->ipo);
    exc = &(st->exc[0]);
    exc2 = &(st->exc2[0]);
    lpi1 = &(st->lpi1);
    lpi2 = &(st->lpi2);
    lpi3 = &(st->lpi3);
    hpi1 = &(st->hpi1);
    hpi2 = &(st->hpi2);
    hpi3 = &(st->hpi3);
    rmso = &(st->rmso_bsynz);

/*                  MAXPIT+MAXORD=166 */
/*  Calculate history scale factor XY and scale filter state */
/* Computing MIN */
    r__1 = *rmso / (*rms + 1e-6f);
    xy = min(r__1,8.f);
    *rmso = *rms;
    i__1 = contrl_1.order;
    for (i__ = 1; i__ <= i__1; ++i__) {
	exc2[i__ - 1] = exc2[*ipo + i__ - 1] * xy;
    }
    *ipo = *ip;
    if (*iv == 0) {
/*  Generate white noise for unvoiced */
	i__1 = *ip;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    exc[contrl_1.order + i__ - 1] = (real) (random_(st) / 64);
	}
/*  Impulse doublet excitation for plosives */
/*       (RANDOM()+32768) is in the range 0 to 2**16-1.  Therefore the
 */
/*       following expression should be evaluated using integers with 
at */
/*       least 32 bits (16 isn't enough), and PX should be in the rang
e */
/*       ORDER+1+0 through ORDER+1+(IP-2) .EQ. ORDER+IP-1. */
	px = (random_(st) + 32768) * (*ip - 1) / 65536 + contrl_1.order + 1;
	r__1 = *ratio / 4 * 1.f;
	pulse = r__1 * 342;
	if (pulse > 2e3f) {
	    pulse = 2e3f;
	}
	exc[px - 1] += pulse;
	exc[px] -= pulse;
/*  Load voiced excitation */
    } else {
	sscale = sqrt((real) (*ip)) / 6.928f;
	i__1 = *ip;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    exc[contrl_1.order + i__ - 1] = 0.f;
	    if (i__ <= 25) {
		exc[contrl_1.order + i__ - 1] = sscale * kexc[i__ - 1];
	    }
	    lpi0 = exc[contrl_1.order + i__ - 1];
	    r__2 = exc[contrl_1.order + i__ - 1] * .125f + *lpi1 * .75f;
	    r__1 = r__2 + *lpi2 * .125f;
	    exc[contrl_1.order + i__ - 1] = r__1 + *lpi3 * 0.f;
	    *lpi3 = *lpi2;
	    *lpi2 = *lpi1;
	    *lpi1 = lpi0;
	}
	i__1 = *ip;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    noise[contrl_1.order + i__ - 1] = random_(st) * 1.f / 64;
	    hpi0 = noise[contrl_1.order + i__ - 1];
	    r__2 = noise[contrl_1.order + i__ - 1] * -.125f + *hpi1 * .25f;
	    r__1 = r__2 + *hpi2 * -.125f;
	    noise[contrl_1.order + i__ - 1] = r__1 + *hpi3 * 0.f;
	    *hpi3 = *hpi2;
	    *hpi2 = *hpi1;
	    *hpi1 = hpi0;
	}
	i__1 = *ip;
	for (i__ = 1; i__ <= i__1; ++i__) {
	    exc[contrl_1.order + i__ - 1] += noise[contrl_1.order + i__ - 1];
	}
    }
/*   Synthesis filters: */
/*    Modify the excitation with all-zero filter  1 + G*SUM */
    xssq = 0.f;
    i__1 = *ip;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k = contrl_1.order + i__;
	sum = 0.f;
	i__2 = contrl_1.order;
	for (j = 1; j <= i__2; ++j) {
	    sum += coef[j] * exc[k - j - 1];
	}
	sum *= *g2pass;
	exc2[k - 1] = sum + exc[k - 1];
    }
/*   Synthesize using the all pole filter  1 / (1 - SUM) */
    i__1 = *ip;
    for (i__ = 1; i__ <= i__1; ++i__) {
	k = contrl_1.order + i__;
	sum = 0.f;
	i__2 = contrl_1.order;
	for (j = 1; j <= i__2; ++j) {
	    sum += coef[j] * exc2[k - j - 1];
	}
	exc2[k - 1] = sum + exc2[k - 1];
	xssq += exc2[k - 1] * exc2[k - 1];
    }
/*  Save filter history for next epoch */
    i__1 = contrl_1.order;
    for (i__ = 1; i__ <= i__1; ++i__) {
	exc[i__ - 1] = exc[*ip + i__ - 1];
	exc2[i__ - 1] = exc2[*ip + i__ - 1];
    }
/*  Apply gain to match RMS */
    r__1 = *rms * *rms;
    ssq = r__1 * *ip;
    gain = sqrt(ssq / xssq);
    i__1 = *ip;
    for (i__ = 1; i__ <= i__1; ++i__) {
	sout[i__] = gain * exc2[contrl_1.order + i__ - 1];
    }
    return 0;
} /* bsynz_ */
@


1.1
log
@Initial revision
@
text
@d3 4
a6 1
$Log$
d11 1
a11 2
extern int bsynz_(real *coef, integer *ip, integer *iv, real *sout, real *rms, real *ratio, real *g2pass);
extern int initbsynz_(void);
d36 4
a39 1
/* $Log: bsynz.f,v $ */
d82 3
a84 2
/* Subroutine */ int bsynz_0_(int n__, real *coef, integer *ip, integer *iv, 
	real *sout, real *rms, real *ratio, real *g2pass)
d88 2
a89 2
    static integer ipo = 0;
    static real rmso = 0.f;
d92 8
a99 28
    static real exc[166] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };
    static real exc2[166] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
	    0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };
    static real lpi0 = 0.f;
    static real lpi2 = 0.f;
    static real lpi3 = 0.f;
    static real hpi0 = 0.f;
    static real hpi2 = 0.f;
    static real hpi3 = 0.f;
d114 1
a114 1
    extern integer random_(void);
d116 1
a116 1
    static real hpi1, lpi1;
d118 4
a121 1
/* $Log: config.fh,v $ */
d135 4
a138 1
/* $Log: contrl.fh,v $ */
d249 10
a258 3
    switch(n__) {
	case 1: goto L_initbsynz;
	}
d263 1
a263 1
    r__1 = rmso / (*rms + 1e-6f);
d265 1
a265 1
    rmso = *rms;
d268 1
a268 1
	exc2[i__ - 1] = exc2[ipo + i__ - 1] * xy;
d270 1
a270 1
    ipo = *ip;
d275 1
a275 1
	    exc[contrl_1.order + i__ - 1] = (real) (random_() / 64);
d285 1
a285 1
	px = (random_() + 32768) * (*ip - 1) / 65536 + contrl_1.order + 1;
d303 6
a308 6
	    r__2 = exc[contrl_1.order + i__ - 1] * .125f + lpi1 * .75f;
	    r__1 = r__2 + lpi2 * .125f;
	    exc[contrl_1.order + i__ - 1] = r__1 + lpi3 * 0.f;
	    lpi3 = lpi2;
	    lpi2 = lpi1;
	    lpi1 = lpi0;
d312 1
a312 1
	    noise[contrl_1.order + i__ - 1] = random_() * 1.f / 64;
d314 6
a319 6
	    r__2 = noise[contrl_1.order + i__ - 1] * -.125f + hpi1 * .25f;
	    r__1 = r__2 + hpi2 * -.125f;
	    noise[contrl_1.order + i__ - 1] = r__1 + hpi3 * 0.f;
	    hpi3 = hpi2;
	    hpi2 = hpi1;
	    hpi1 = hpi0;
a365 31
/*   Print test data */
/* 	IF(LISTL.GE.5) THEN */

/*           I changed the range of indices to print within NOISE from */
/*           1,IP+ORDER to 1+ORDER,IP+ORDER since indices 1 through ORDER 
*/
/*           of NOISE are never used.  This avoids printing out their */
/*           "garbage" values. */

/* 	   IF(IV.NE.0) */
/*     1     WRITE(FDEBUG,980) 'NOISE:',(NOISE(I),I=1+ORDER,IP+ORDER) */
/* 	   WRITE(FDEBUG,980) 'EXC:',  (EXC(I),  I=1,IP+ORDER) */
/* 	   WRITE(FDEBUG,980) 'EXC2:', (EXC2(I), I=1,IP+ORDER) */
/* 	   WRITE(FDEBUG,980) 'SOUT:', (SOUT(I), I=1,IP) */
/* 980	   FORMAT(1X,A,100(/1X,10F10.1)) */
/* 	END IF */
    return 0;

L_initbsynz:
    ipo = 0;
    for (i__ = 1; i__ <= 166; ++i__) {
	exc[i__ - 1] = 0.f;
	exc2[i__ - 1] = 0.f;
    }
    lpi0 = 0.f;
    lpi2 = 0.f;
    lpi3 = 0.f;
    hpi0 = 0.f;
    hpi2 = 0.f;
    hpi3 = 0.f;
    rmso = 0.f;
a367 13

/* Subroutine */ int bsynz_(real *coef, integer *ip, integer *iv, real *sout, 
	real *rms, real *ratio, real *g2pass)
{
    return bsynz_0_(0, coef, ip, iv, sout, rms, ratio, g2pass);
    }

/* Subroutine */ int initbsynz_(void)
{
    return bsynz_0_(1, (real *)0, (integer *)0, (integer *)0, (real *)0, (
	    real *)0, (real *)0, (real *)0);
    }

@