/**********************************************************************
MPEG-4 Audio VM Module
parameter based codec - HILN: harmonic/individual lines and noise
                              extraction



This software module was originally developed by

Heiko Purnhagen (University of Hannover / Deutsche Telekom Berkom)
Bernd Edler (University of Hannover / Deutsche Telekom Berkom)

and edited by

in the course of development of the MPEG-2 NBC/MPEG-4 Audio standard
ISO/IEC 13818-7, 14496-1,2 and 3. This software module is an
implementation of a part of one or more MPEG-2 NBC/MPEG-4 Audio tools
as specified by the MPEG-2 NBC/MPEG-4 Audio standard. ISO/IEC gives
users of the MPEG-2 NBC/MPEG-4 Audio standards free license to this
software module or modifications thereof for use in hardware or
software products claiming conformance to the MPEG-2 NBC/ MPEG-4 Audio
standards. Those intending to use this software module in hardware or
software products are advised that this use may infringe existing
patents. The original developer of this software module and his/her
company, the subsequent editors and their companies, and ISO/IEC have
no liability for use of this software module or modifications thereof
in an implementation. Copyright is not released for non MPEG-2
NBC/MPEG-4 Audio conforming products. The original developer retains
full right to use the code for his/her own purpose, assign or donate
the code to a third party and to inhibit third party from using the
code for non MPEG-2 NBC/MPEG-4 Audio conforming products. This
copyright notice must be included in all copies or derivative works.

Copyright (c) 1997.


NOTE:
=====
This module provides the basic funtionality to generate all elements of
a bit stream for HILN parameter based coding. However it is not fully
optimised for providing the best possible audio quality.


Source file: indilinextr.c

$Id: indilinextr_basic.c,v 1.6 1998/04/30 18:30:21 purnhage Exp $


Authors:
HP    Heiko Purnhagen, Uni Hannover <purnhage@tnt.uni-hannover.de>
BE    Bernd Edler, Uni Hannover <edler@tnt.uni-hannover.de>

Changes:
13-nov-97   HP/BE first version based on advanced indilinextr.c
27-nov-97   HP    fixes for some C compilers (no empty struct's)
16-feb-98   HP    fixed longwin seg fault bug (thanks to Matti Koskinen)
16-feb-98   HP    improved noiseFreq, fixed harm maxLineFreq
01-mar-98   HP    ILXharmCorrEnable
31-mar-98   HP    completed IndiLineExtractFree()
09-apr-98   HP    ILXconfig
29-apr-98   HP    fixed specPwr calc
**********************************************************************/


#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include "libuhd_psy.h"		/* psychoacoustic model */
#include "indilinecom.h"	/* indiline common module */
#include "indilinextr.h"	/* indiline extraction module */

#include "uhd_fft.h"		/* FFT and ODFT module */

/* units & formats: */
/*  Freq: frequency in Hz */
/*  Ampl: peak amplitude of sine wave */
/*  Phase: phase in rad at frame center */
/*  Index: NumBit LSBs, to be transmitted MSB first */
/*  Pred: line predecessor index: 0=no pred., 1=line[0], 2=line[1], ... */


/* ---------- declarations (global) ---------- */

#define min(a,b) ((a) < (b) ? (a) : (b))
#define max(a,b) ((a) > (b) ? (a) : (b))

#define IL_MINFLOAT ((float)1e-35)
#define ATAN2M_MIN_P 1e-6


/* debug flags */

#define MINENVAMPL 0.1

/* window functions */
#define SHORT_RAT 8			/* ratio for short overlap */

#define LONG_THR 5

#define FMAX_FSAMPLE_RATIO 0.49


/* ---------- declarations (internal, for IndiLineExtract) ---------- */

/* parameters for array declarations */

/* line continuation and line parameter limits */
#define DFMAX	(.05)		/* max. freq. change for continued lines */
#define DFQMAX	(.15)		/* max. quant. freq. change for cont. lines */
#define FAMAX	(4.)		/* max. ampl. factor for continued lines */
#define FAMIN	(1./FAMAX)	/* min. ampl. factor for continued lines */

#define ILX_FMIN (30.0)		/* min. frequency of extracted line */
#define ILX_MINAMPL (0.1/32768)	/* min. amplitude of extracted line */
				/* 1.0/32768 doesn't give enough lines */
#define ILX_HARMMINAMPL (0.01)	/* count harm lines above -40dB re. max harm */
#define ILX_MINNUMHARM (3)
#define MAX_HARMAMPLPEAK (5.)

#define MAXHARMINDI (3)
#define HARMSPECTHRES (0.5)		/* 0.7 */
#define HARMPWRTHRES (0.8)
#define HARMCORRTHRES (0.9)
#define LINESPECTHRES (0.5)		/* 0.7 */

/* ---------- declarations for envelope and window switching ---------- */

#define ENV_OVER_F (1./16.)	/* overlap factor to find as,ae */
				/* should be 1/(2*OV_RAT) in asdec_o */
#define SINGLE_RATE_THR (1)	/* rate threshold for single slope envelope */

/* parameters for hilbert transform and filter */
#define HWINLEN2_F (1./16.)	/* hilbert window length */
				/* (-HWINLEN2_F*FRMLEN..HWINLEN2_F*FRMLEN) */
#define HWINHAN (1)		/* 0=rectangular 1=hanning window */
#define FWINLEN2_F (1./16.)	/* filter window length (0 -> filter off) */
				/* (-FWINLEN2_F*FRMLEN..FWINLEN2_F*FRMLEN) */
#define FWINHAN (1)		/* 0=rectangular 1=hanning window */

#define ENV_MINFREQ (100.)

/* ---------- declarations for harmonic component extraction ---------- */

#define SPECZEROPAD (2)		/* spectrum zero padding factor (1=off) */
#define SPECMAXWIN (4000)	/* max spectrum window width [Hz] */
#define SPECMINPWR (1e-10)	/* min power spectrum value */
#define MAXFUNDFREQ (1500)	/* max fundamental frequency [Hz] */

static float ILXHminFundFreq = 20.0;
static float ILXHmaxFundFreq = 20000.0;


/* ---------- typedef's for data structures ---------- */

typedef struct ILXEstatusStruct ILXEstatus;	/* ILXE status handle */

typedef struct ILXHstatusStruct ILXHstatus;	/* ILXH status handle */

typedef struct ILXLstatusStruct ILXLstatus;	/* ILXL status handle */

/* ---------- declarations (data structures) ---------- */

/* status variables and arrays for envelope parameter calculation */

struct ILXEstatusStruct {
  /* general parameters */
  int debugLevel;		/* debug level (0=off) */
  int frameLen;			/* frame length */
  int envOver;			/* overlap for envelope calculation */
  int hwinlen2;			/* half of hilbert transform filter length */
  int fwinlen2;			/* half of band pass filter length */
  int envLen;			/* size of env array */
  /* pointers to fixed content arrays ("static") */
  float *hilb;			/* impulse response hilbert transform */
  float *filt;			/* impulse response band pass filter */
  /* pointers to variable content arrays */
  float *xfilt;			/* filter output */
  float *env;			/* output of filter and hilbert transform */
};

/* status variables and arrays for harmonic tone parameter calculation */

struct ILXHstatusStruct {
  /* general parameters */
  int debugLevel;		/* debug level (0=off) */
  int winLen;			/* spectrum window length */
  int specLen;			/* total spectrum length */
  float specWidth;		/* total spectrum bandwidth */

  /* fixed content arrays */
  float *win;			/* frequency window */

  /* variable content arrays */
  float *re;			/* cepstrum (real) */
  float *im;			/* cepstrum (imag.) */
};

/* status variables and arrays for line parameter calculation */

struct ILXLstatusStruct {
  int dummy;			/* some compilers don't like empty struct's */
};

/* status variables and arrays for individual line extraction (ILX) */

struct ILXstatusStruct		/* ILX status handle struct */
{
  /* general parameters */
  int debugLevel;		/* debug level (0=off) */
  int frameLen;			/* frame length */
  float fSample;		/* sampling frequency */
  float maxLineFreq;		/* max. allowed frequency */
  float maxLineAmpl;		/* max. allowed amplitude */
  int maxNumLine;		/* maximum nuber of lines */
  int maxNumEnv;
  int maxNumHarm;
  int maxNumHarmLine;
  int maxNumNoisePara;
  float maxNoiseFreq;
  int maxNumAllLine;
  int bufLen;			/* buffer length */

  /* parameters for frequency and amplitude estimation and psychoacoustics */
  int winLen;			/* window length  */
  int fftLen;			/* transform length  */
  int addStart;			/* add. samples before center window */

  /* fixed content arrays */
  float *costab;		/* cosine table for dft coefficients */
  float *longWin;		/* long window (full overlap) */

  /* arrays for return values */
  float **envPara;		/* envelope parameters */
  float *lineFreq;		/* line frequencies [Hz] */
  float *lineAmpl;		/* line amplitudes */
  float *linePhase;		/* line phase values */
  int *lineEnv;			/* line envelope flags */
  int *linePred;		/* line predecessor indices */
  int *numHarmLine;
  float *harmFreq;
  float *harmFreqStretch;
  int *harmLineIdx;
  int *harmEnv;
  int *harmPred;
  float *harmRate;
  float *noisePara;

  /* variable content arrays */
  float *cos;			/* cosine with current constant frequency */
  float *sin;			/* sine with current constant frequency */
  float *envWin;		/* window with current envelope */
  float *x;			/* residual signal */
  float *xOrgRe;		/* real in/out of dft */
  float *xOrgIm;		/* imag in/out of dft */
  float *xmodRe;		/* real modulated residual */
  float *xmodIm;		/* imag modulated residual */
  float *xsyn;			/* accumulated synth. signal (all lines) */
  float *xsynRe;		/* spectrum of synth. signal (real part) */
  float *xsynIm;		/* spectrum of synth. signal (imag part) */
  float *lineDf;		/* freq. change within the frame */
  int *lflag;			/* flag for excluding lines from freq. est. */
  int *prevLineCont;		/* index of prev. lines for continuation */
  float *harmCosFact;
  float *harmSinFact;
  float *harmFc;
  float *harmFe;
  int *harmPsy;
  
  /* PSYCHOACOUSTIC, currently under test */
  float *excitSynthFft; /* FFT-related excitation of synth. signal */
 
  /* frame-to-frame memory */
  int prevNumEnv;
  float **prevEnvPara;		/* envelope parameters */
  int prevNumLine;		/* number of extracted lines */
  float *prevLineFreq;		/* line frequencies [Hz] */
  float *prevLineAmpl;		/* line amplitudes */
  float *prevLinePhase;		/* line phase values */
  int *prevLineEnv;		/* line envelope flags */
  int prevNumHarm;
  int *prevNumHarmLine;
  float *prevHarmFreq;
  float *prevHarmFreqStretch;
  int *prevHarmLineIdx;
  int *prevHarmEnv;

  float prevFundaF;

  /* status handles */
  ILXEstatus *envStatus;		/* envelope parameter estimation */
  ILXLstatus *lineStatus;		/* line parameter estimation */
  ILXHstatus *harmStatus;		/* harm. tone parameter estimation */
};


/* ---------- functions (internal) ---------- */


/* atan2m  */

static double atan2m(double im, double re)
{
  if(fabs(im)<ATAN2M_MIN_P && fabs(re)<ATAN2M_MIN_P)
    return(0);
  else if(fabs(re)<ATAN2M_MIN_P && im>0)
    return(M_PI/2);
  else if(fabs(re)<ATAN2M_MIN_P && im<0)
    return(-M_PI/2);
  else
    return(atan2(im,re));
}


/* ILXEnvInit() */
/* Init envelope parameter estimation. */

static ILXEstatus *ILXEnvInit (
  int frameLen,			/* in: samples per frame */
  int *addSampStart,		/* out: add. samp. at start */
  int *addSampEnd,		/* out: add. samp. at end */
  int debugLevel)		/* in: debug level (0=off) */
				/* returns: ILXE status handle */
{
  ILXEstatus	*ILXE;
  int i,envAdd,xfiltLen;
  float temp;

  if ((ILXE = (ILXEstatus *) malloc(sizeof(ILXEstatus)))==NULL)
    IndiLineExit("ILXEnvInit: memory allocation error");

  ILXE->debugLevel = debugLevel;
  ILXE->frameLen = frameLen;
  ILXE->envOver = (int)(ENV_OVER_F*frameLen);
  ILXE->hwinlen2 = (int)(HWINLEN2_F*frameLen);
  ILXE->fwinlen2 = (int)(FWINLEN2_F*frameLen);
  envAdd = max(ILXE->envOver,ILXE->frameLen/2);	/* add. samples start & end */
  ILXE->envLen = frameLen+2*envAdd;
  *addSampStart = *addSampEnd = envAdd+ILXE->hwinlen2+ILXE->fwinlen2;

  /* fixed content arrays */
  if(((ILXE->hilb = (float *)malloc((2*ILXE->hwinlen2+1)*sizeof(float)))==NULL)
  || ((ILXE->filt = (float *)malloc((2*ILXE->fwinlen2+1)*sizeof(float)))==NULL))
    IndiLineExit("ILXEnvInit: memory allocation error");

  /* variable content arrays */
  xfiltLen = ILXE->envLen+2*ILXE->hwinlen2;
  if(((ILXE->xfilt = (float *)malloc(xfiltLen*sizeof(float)))==NULL)
  || ((ILXE->env = (float *)malloc(ILXE->envLen*sizeof(float)))==NULL))
    IndiLineExit("ILXEnvInit: memory allocation error");

  /* init hilbert and filter */
  ILXE->hilb[ILXE->hwinlen2] = 0;
  for (i=1;i<=ILXE->hwinlen2;i++) {
    if (i%2) {
      ILXE->hilb[ILXE->hwinlen2+i] = 1/(M_PI/2*i);
      if (HWINHAN)
	ILXE->hilb[ILXE->hwinlen2+i] *= .5+.5*cos(M_PI/ILXE->hwinlen2*i);
    }
    else
      ILXE->hilb[ILXE->hwinlen2+i] = 0;
    ILXE->hilb[ILXE->hwinlen2-i] = -ILXE->hilb[ILXE->hwinlen2+i];
  }
  if (ILXE->fwinlen2==0)
    ILXE->filt[0] = 1;
  else {
    ILXE->filt[ILXE->fwinlen2] = (ILXE->frameLen/2-1)/(float)(ILXE->frameLen/2);
    temp = 0;
    for (i=1;i<=ILXE->fwinlen2;i++) {
      if (i%2==0) {
	ILXE->filt[ILXE->fwinlen2+i] = -1/(float)(ILXE->frameLen/2);
	if (FWINHAN)
	  ILXE->filt[ILXE->fwinlen2+i] *= .5+.5*cos(M_PI/ILXE->fwinlen2*i);
      }
      else
	ILXE->filt[ILXE->fwinlen2+i] = 0;
      ILXE->filt[ILXE->fwinlen2-i] = ILXE->filt[ILXE->fwinlen2+i];
      temp += 2*ILXE->filt[ILXE->fwinlen2+i];
    }
    /* make FILT(f=0) -> 0 */
    temp = -ILXE->filt[ILXE->fwinlen2]/temp;
    for (i=1;i<=ILXE->fwinlen2;i++) {
      ILXE->filt[ILXE->fwinlen2+i] *= temp;
      ILXE->filt[ILXE->fwinlen2-i] *= temp;
    }
    /* make FILT(f=1/4) -> 1 */
    temp = ILXE->filt[ILXE->fwinlen2];
    for (i=2;i<=ILXE->fwinlen2;i+=2) {
      if (i%4==0)
	temp += 2*ILXE->filt[ILXE->fwinlen2+i];
      else
	temp -= 2*ILXE->filt[ILXE->fwinlen2+i];
    }
    for(i=0;i<ILXE->fwinlen2*2+1;i++)
      ILXE->filt[i] /= temp;
  }

  return(ILXE);
}


/* ILXEnvFree() */
/* Free memory allocated by ILXEnvInit(). */

static void ILXEnvFree (
  ILXEstatus *ILXE)		/* in: ILXE status handle*/
{
  free(ILXE->hilb);
  free(ILXE->filt);
  free(ILXE->xfilt);
  free(ILXE->env);
  free(ILXE);
}


/* ILXLineInit() */
/* Init line parameter estimation. */

static ILXLstatus *ILXLineInit (
  int debugLevel)		/* in: debug level (0=off) */
				/* returns: ILXL status handle */
{
  return NULL;
}


/* ILXLineFree() */
/* Free memory allocated by ILXLineInit(). */

static void ILXLineFree (
  ILXLstatus *ILXL)		/* in: ILXL status handle*/
{
}


/* ILXHarmInit() */
/* Init harmonic tone parameter estimation. */

static ILXHstatus *ILXHarmInit (
  int debugLevel,		/* in: debug level (0=off) */
  int specLen,			/* in: spectrum length */
  float specWidth)		/* in: spectrum bandwidth [Hz] */
				/* returns: ILXH status handle */
{
  ILXHstatus	*ILXH;
  int i,len;

  if ((ILXH = (ILXHstatus *) malloc(sizeof(ILXHstatus)))==NULL)
    IndiLineExit("ILXHarmInit: memory allocation error");

  ILXH->debugLevel = debugLevel;
  ILXH->winLen = (int)(min(specWidth,SPECMAXWIN)/specWidth*specLen+.5);
  len = 1;
  while (len < ILXH->winLen)
    len *= 2;
  ILXH->specLen = len*SPECZEROPAD*2;
  ILXH->specWidth = (specWidth/specLen)*ILXH->specLen;

  if(((ILXH->win = (float *)malloc(ILXH->winLen*sizeof(float)))==NULL)
  || ((ILXH->re = (float *)malloc(ILXH->specLen*sizeof(float)))==NULL)
  || ((ILXH->im = (float *)malloc(ILXH->specLen*sizeof(float)))==NULL))
    IndiLineExit("ILXHarmInit: memory allocation error");

  for (i=0; i<ILXH->winLen; i++)
    ILXH->win[i] = (1+cos((i+.5)/(ILXH->winLen)*M_PI))/2;

  return ILXH;
}


/* ILXHarmFree() */
/* Free memory allocated by ILXHarmInit(). */

static void ILXHarmFree (
  ILXHstatus *ILXH)		/* in: ILXH status handle*/
{
  free(ILXH->win);
  free(ILXH->im);
  free(ILXH->re);
  free(ILXH);
}


/* ILXEnvEstim() */
/* Estimate envelope parameter. */

static int ILXEnvEstim (
  ILXEstatus *ILXE,		/* in: ILXE status handle */
  float *x,			/* in: input signal */
  float *envPara,		/* out: envelope parameters */
  float *envWin,		/* out: env. built from par. */
  int envWinLen)		/* length of env. window */
				/* returns: envFlag */
{
  int i,ii,im,maxi,envFlag,envAdd,envWinOff;
  float rms,am,r0,r1,r2,rx0,rx1,ra,rb,wght,temp;
  float atk_ratei,dec_ratei;
  float *xfilt,*filt,*hilb,*env;

  /* set pointers */
  envAdd = max(ILXE->envOver,ILXE->frameLen/2);	/* add. samples start & end */
  xfilt = &ILXE->xfilt[envAdd+ILXE->hwinlen2];	/* filter output */
  filt = &ILXE->filt[ILXE->fwinlen2];		/* filter coefficients */
  hilb = &ILXE->hilb[ILXE->hwinlen2];		/* hilbert coefficients */
  env = &ILXE->env[envAdd];			/* output of filt. and hilb. */

  /* approximate signal envelope based on filter and hilbert transform */
  for (i=-envAdd-ILXE->hwinlen2;i<ILXE->frameLen+envAdd+ILXE->hwinlen2;i++){
    temp = 0;
    for (ii=-ILXE->fwinlen2;ii<=ILXE->fwinlen2;ii++)
      temp += filt[ii]*x[i-ii];
    xfilt[i] = temp;
  }
  for (i=-envAdd;i<ILXE->frameLen+envAdd;i++){
    temp = 0;
    for (ii=-ILXE->hwinlen2;ii<=ILXE->hwinlen2;ii++)
      temp += hilb[ii]*xfilt[i-ii];
    env[i] = sqrt(xfilt[i]*xfilt[i]+temp*temp);
  }

  /* calc max and rms */
  rms = 0;
  am = 1;
  im = 0;
  for(i=-ILXE->envOver;i<ILXE->frameLen+ILXE->envOver;i++) {
    rms += env[i]*env[i];
    if (env[i]>am) {
      am = env[i];
      im = i;
    }
  }
  im = max(0,min(ILXE->frameLen-1,im));
  rms = sqrt(rms/(ILXE->frameLen+2*ILXE->envOver))/am;

  /* calc envelope parameters using regression */

  /* attack */
  r2 = rx1 = 0;
  for(i=0;i<max(im+ILXE->envOver,ILXE->frameLen/2);i++) {
    temp = env[im-i]/am;
    if (temp<=rms)
      wght = 0;
    else
      wght = pow((temp-rms)/(1-rms),1+i/(float)(ILXE->frameLen/4))*
	(1+i/(float)(ILXE->frameLen/4));
    /* wght = pow((temp-rms)/(1-rms),
       max(1/3.,min(3,i/(float)(ILXE->frameLen/2))))*
       (1+i/(float)(ILXE->frameLen/2)); */
    temp = 1-temp;
    r2 += wght*i*i;
    rx1 += wght*temp*i;
  }
  /* line env[i]=1-ra*(im-i) */
  if (r2==0)
    ra = 0;
  else
    ra = rx1/r2;
  if (ILXE->debugLevel)
    printf("atk: ra=%f",ra);
  atk_ratei = max(0,ra);
  /* calc tmax */;
  maxi = im;

  /* decay */
  r2 = rx1 = 0;
  for(i=0;i<max(ILXE->frameLen+ILXE->envOver-maxi,ILXE->frameLen/2);i++) {
    temp = env[maxi+i]/am;
    if (temp<=rms)
      wght = 0;
    else
      wght = pow((temp-rms)/(1-rms),1+i/(float)(ILXE->frameLen/4))*
	(1+i/(float)(ILXE->frameLen/4));
    /* wght = pow((temp-rms)/(1-rms),
       max(1/3.,min(3,i/(float)(ILXE->frameLen/2))))*
       (1+i/(float)(ILXE->frameLen/2)); */
    temp = 1-temp;
    r2 += wght*i*i;
    rx1 += wght*temp*i;
  }
  /* line env[i]=1-ra*(i-maxi) */
  if (r2==0)
    ra = 0;
  else
    ra = rx1/r2;
  if (ILXE->debugLevel)
    printf("   dec: ra=%f\n",ra);
  dec_ratei = max(0,ra);

  if (ILXE->debugLevel) {
    printf("rms=%f am=%f im=%d\n",rms,am,im);
    printf("prelim. atk_ratei=%f dec_ratei=%f maxi=%d (%f)\n",
	   atk_ratei,dec_ratei,maxi,maxi/(float)ILXE->frameLen);
  }
  /* check if single slope envelope */
  if ((atk_ratei*ILXE->frameLen < SINGLE_RATE_THR) &&
      (dec_ratei*ILXE->frameLen < SINGLE_RATE_THR)) {
    r0 = r1 = r2 = rx0 = rx1 = 0;
    for(i=-ILXE->envOver;i<=ILXE->frameLen+ILXE->envOver;i++) {
      temp = env[i]/am;
      if (temp<rms)
	wght = 0;
      else
	wght = (temp-rms)/(1-rms);
      r0 += wght;
      r1 += wght*i;
      r2 += wght*i*i;
      rx0 += wght*temp;
      rx1 += wght*temp*i;
    }
    /* line: env[i]=ra*i+rb */
    temp = (r1*r1-r2*r0);
    if (temp==0)
      ra = rb = 0;
    else {
      ra = (rx0*r1-rx1*r0)/temp;
      rb = (r1*rx1-r2*rx0)/temp;
    }
    if (ILXE->debugLevel)
      printf("single slope: ra=%f rb=%f",ra,rb);
    if (fabs(ra*ILXE->frameLen) < 0.5) {
      atk_ratei = dec_ratei = 0;
      maxi = 0;
      if (ILXE->debugLevel)
	printf("no env\n");
    }
    else {
      if (ra >= 0) {
	atk_ratei = ra;
	dec_ratei = 0;
	maxi = ILXE->frameLen-1;
      }
      else {
	atk_ratei = 0;
	dec_ratei = -ra;
	maxi = 0;
      }
      if (ILXE->debugLevel)
	printf(" (maxi=%d)\n",maxi);
    }
  }

  envFlag = (atk_ratei != 0) || (dec_ratei != 0);

  if (ILXE->debugLevel)
    printf("env=%d atk_ratei=%f dec_ratei=%f maxi=%d (%f)\n",
	   envFlag,atk_ratei,dec_ratei,maxi,maxi/(float)ILXE->frameLen);

  /* build window according to envelope */

  envWinOff = (envWinLen-ILXE->frameLen)/2;
  for(i=0;i<envWinLen;i++)
    envWin[i] = 0;

  i = maxi+envWinOff;
  while (i>=0) {
    envWin[i] = 1-(atk_ratei)*(maxi-(i-envWinOff));
    if (envWin[i] <= 0) {
      envWin[i] = 0;
      i = -1;
    } 
    else
      i--;
  }
  i = maxi+envWinOff;
  while (i<envWinLen) {
    envWin[i] = 1-(dec_ratei)*((i-envWinOff)-maxi);
    if (envWin[i] <= 0) {
      envWin[i] = 0;
      i = envWinLen;
    } 
    else
      i++;
  }

  envPara[0] = maxi/(float)ILXE->frameLen;
  envPara[1] = atk_ratei*(float)ILXE->frameLen;
  envPara[2] = dec_ratei*(float)ILXE->frameLen;

  return(envFlag);
}


/* FindPeak() */
/* Find peak (local maximum) in array. */
/* x[p-2] < x[p-1] < x[p] > x[p+1] > x[p+2] */

static int FindPeak (
  float *x,			/* in: array */
  float xmin,			/* in: min value of peak (treshold) */
  int imin,			/* in: min index (imin+2<=p) */
  int imax)			/* in: max index (p<=imax-2) */
				/* returns: index p (or -1 if not found) */
{
  int i,ipeak;
  int state;

  ipeak = -1;
  state = 0;
  for (i=imin; i<imax; i++) {
    if (x[i] < x[i+1])
      switch (state) {
	case 0: state=1; break;
	case 1: state=2; break;
	case 2: state=2; break;
	case 3: state=0; break;
      }
    else
      switch (state) {
	case 0: state=0; break;
	case 1: state=0; break;
	case 2: state=3; break;
	case 3:
	  state=0;
	  if (x[i-1] > xmin) {
	    ipeak = i-1;
	    xmin = x[i-1];
	  }
	  break;
      }
  }
  return ipeak;
}


/* ILXHarmFundFreq() */
/* Estimate fundamental frequency (maximum in real cepstrum). */

static float ILXHarmFundFreq (
  ILXHstatus *ILXH,		/* in: ILXH status handle*/
  float *xRe,			/* in: ODFT spectrum (real) */
  float *xIm)			/* in: ODFT spectrum (imag.) */
				/* returns: fundamental frequency [Hz] */
{
  int i,i0;
  float idx,idxmax;
  float cepmax;
  int t,tt;
  int i1,i2;

  if (ILXHminFundFreq == ILXHmaxFundFreq)
    return ILXHmaxFundFreq;		/* no funda freq estim. */

  for (i=0; i<ILXH->winLen; i++) {
    ILXH->re[i] = ILXH->re[ILXH->specLen-1-i] =
      log(max(SPECMINPWR,xRe[i]*xRe[i]+xIm[i]*xIm[i]))*ILXH->win[i];
    ILXH->im[i] = ILXH->im[ILXH->specLen-1-i] = 0;
  }
  for (i=ILXH->winLen; i<ILXH->specLen/2; i++) {
    ILXH->re[i] = ILXH->re[ILXH->specLen-1-i] = 0;
    ILXH->im[i] = ILXH->im[ILXH->specLen-1-i] = 0;
  }
  UHD_iodft(ILXH->re,ILXH->im,ILXH->specLen);

  i0 = (int)(ILXH->specWidth/MAXFUNDFREQ+.5);
  i1 = (int)(ILXH->specWidth/ILXHmaxFundFreq+.5);
  i2 = (int)(ILXH->specWidth/ILXHminFundFreq+.5);

  if (ILXH->debugLevel) {
    printf("0.123456789 ");
    for (i=0; i<ILXH->specLen/2; i++)
      printf("%7f ",ILXH->re[i]);
  }

  i = FindPeak(ILXH->re,0,max(i0-2,i1-2),min(ILXH->specLen/2,i2+2));
  cepmax = ILXH->re[i];
  idxmax = i;
  t = 1;
  for (tt=2; tt<=5; tt++) {
    i = (int)(idxmax*tt);
    if (i+3 < ILXH->specLen/2 &&
	(i = FindPeak(ILXH->re,cepmax/2,i-3,i+3)) != -1) {
      idxmax = (float)i/tt;
      t = tt;
    }
  }
  for (tt=4; tt>1; tt--) {
    idx = idxmax/tt;
    i = (int)(idx+.5);
    if (i >= i1 &&
	FindPeak(ILXH->re,cepmax/2,i-3,i+3) != -1) {
      idxmax = idx;
      t = -tt;
    }
  }

  if (ILXH->debugLevel) {
    printf("%7f %d\n",idxmax,t);
    printf("fndfrq=%7.1f  cep0=%7f  cepmax=%7f  idxmax=%5.1f  r=1/%1d\n",
	   ILXH->specWidth/idxmax,*ILXH->re,cepmax,idxmax,t);
    printf("98765 %7.1f %7f %7f %5.1f %1d\n",
	   ILXH->specWidth/idxmax,*ILXH->re,cepmax,idxmax,t);
  }

  return ILXH->specWidth/idxmax;
}


/* CalcFactors */
/* calculate factors for spectral line with and without envelope, */
/* calc. reduction of signal variance */

static void CalcFactors(
	float	*x,			/* in: input signal */
	float	*cos,			/* in: cosine for given line */
	float	*sin,			/* in: sine for given line */
	float	*win,			/* in: window for corr. and error */
	float	*envWin,		/* in: add. env. (NULL==no env.) */
	int	winLen,			/* in: window length */
	float	*cosFact,		/* out: factor for cosine */
	float	*sinFact,		/* out: factor for sine */
	float	*qeDiff,		/* out: neg. signal power */
	float	*cosFactEnv,		/* out: factor for cosine (env) */
	float	*sinFactEnv,		/* out: factor for sine (env) */
	float	*qeDiffEnv,		/* out: neg. signal power (env) */
					/* use env if qeDiffEnv<qeDiff */
	int	debugLevel)		/* in: debug level (0=off) */
{
  int i;
  float xCos,xSin,wi_2,wi_ex;
  float ccSum,ssSum,csSum,cxSum,sxSum;
  float ccSumEnv,ssSumEnv,csSumEnv,cxSumEnv,sxSumEnv;
  float det;

  ccSum = ssSum = csSum = cxSum = sxSum = 0;
  ccSumEnv = ssSumEnv = csSumEnv = cxSumEnv = sxSumEnv = 0;
  for(i=0;i<winLen;i++) {
    xCos = x[i]*cos[i];
    xSin = x[i]*sin[i];
    wi_2 = win[i]*win[i];
    ccSum += wi_2*cos[i]*cos[i];
    ssSum += wi_2*sin[i]*sin[i];
    csSum += wi_2*cos[i]*sin[i];
    cxSum += wi_2*xCos;
    sxSum += wi_2*xSin;
    if(envWin!=NULL) {
      wi_2 = (envWin[i]*win[i])*(envWin[i]*win[i]);
      wi_ex = (envWin[i]*win[i])*win[i];
      ccSumEnv += wi_2*cos[i]*cos[i];
      ssSumEnv += wi_2*sin[i]*sin[i];
      csSumEnv += wi_2*cos[i]*sin[i];
      cxSumEnv += wi_ex*xCos;
      sxSumEnv += wi_ex*xSin;
    }
  }
  det = max(ccSum*ssSum-csSum*csSum,IL_MINFLOAT);
  *cosFact = (ssSum*cxSum-csSum*sxSum)/det;
  *sinFact = (-ccSum*sxSum+csSum*cxSum)/det;
  *qeDiff = -(*cosFact)*cxSum+(*sinFact)*sxSum;

  if (debugLevel)
    printf("noenv: a=%7.1f, de=%7.1e",
	   sqrt((*cosFact)*(*cosFact)+(*sinFact)*(*sinFact)),*qeDiff);

  if(envWin!=NULL) {
    det = max(ccSumEnv*ssSumEnv-csSumEnv*csSumEnv,IL_MINFLOAT);
    *cosFactEnv =(ssSumEnv*cxSumEnv-csSumEnv*sxSumEnv)/det;
    *sinFactEnv =(-ccSumEnv*sxSumEnv+csSumEnv*cxSumEnv)/det;
    *qeDiffEnv = -(*cosFactEnv)*cxSumEnv+(*sinFactEnv)*sxSumEnv;

    if (debugLevel)
      printf(" env: a=%7.1f, de=%7.1e",
	     sqrt((*cosFactEnv)*(*cosFactEnv)+(*sinFactEnv)*(*sinFactEnv)),
	     *qeDiffEnv);
  }

  if (debugLevel)
    printf("\n");
}


/* WinOdft */
/* multiply with window function and calc. ODFT */

static void WinOdft (
  float	*x,			/* in: input signal */
  float	*win,			/* in: window function */
  int	winLen,			/* in: length of input vector and window */
  int	fftLen,			/* in: transform length */
  float	*xReal,			/* out: transformed vector (real part) */
  float	*xImag)			/* out: transformed vector (imag. part) */
{
  int i;

  if(winLen>fftLen)
    IndiLineExit("WinOdft: winLen>fftLen");

  for(i=0;i<winLen;i++) {
    xReal[i] = win[i]*x[i];
    xImag[i] = 0;
  }
  for(i=winLen;i<fftLen;i++)
    xReal[i] = xImag[i] = 0;
  UHD_odft(xReal,xImag,fftLen);
}


/* EnvPwrFact() */
/* calculate power ratio for envelope */

static float EnvPwrFact (
  float *envPara)
{
  float tAtk, tDec;
  float pAtk, pDec, ac1, ac2;
    
  tAtk = min(envPara[0], 1./max(1.,envPara[1]));
  tDec = min(1.-envPara[0], 1./max(1.,envPara[2]));
  
  ac1 = 1.-envPara[1]*tAtk;
  ac2 = 1.-envPara[2]*tDec;
  
  pAtk = (envPara[1]==0) ? tAtk : 1./(3.*envPara[1])*(1.-ac1*ac1*ac1);
  pDec = (envPara[2]==0) ? tDec : 1./(3.*envPara[2])*(1.-ac2*ac2*ac2);

  return(pAtk+pDec);
}


/* NoiseDct */
/* DCT of arbitrary length (sizex), calculates */
/* only the sizey first output values */

static void NoiseDct(
  float *x,			/* in:  DCT input */
  float *y,			/* out: DCT output */
  int sizex,			/* in:  length of vector x */
  int sizey,			/* in:  length of vector y */
  int sizedct)			/* in:  length of DCT */
{
  int i,j;

  for(i=0;i<sizey;i++) {
    y[i] = 0;
    for(j=0;j<sizex;j++)
      y[i] += x[j]*cos(M_PI*(float)i*((float)j+.5)/(float)sizedct);
    y[i] *= sqrt(2./(float)sizedct);
  }
  y[0] /= sqrt(2.);
}


/* ---------- functions (global) ---------- */

/* IndiLineExtractInit */
/* init for IndiLineExtract functions */

ILXstatus *IndiLineExtractInit (
  int frameLen,			/* in: samples per frame */
  float fSample,		/* in: sampling frequency */
  float maxLineFreq,		/* in: max line frequency */
  float maxLineAmpl,		/* in: max line amplitude */
  int maxNumLine,		/* in: max num lines */
  int maxNumEnv,		/* in: max num envelopes */
  int maxNumHarm,		/* in: max num harmonic tones */
  int maxNumHarmLine,		/* in: max num lines of harmonic tone */
  int maxNumNoisePara,		/* in: max num noise parameters */
  float maxNoiseFreq,		/* in: max noise frequency */
  ILXconfig *cfg,		/* in: ILX configuration */
				/*     or NULL */
  int debugLevel)		/* in: debug level (0=off) */
				/* returns: ILX status handle */
{
  ILXstatus *ILX;
  int i,imax_psy,winLen;
  int addStartEnv,addEndEnv;
  float pf;

  if ((ILX = (ILXstatus *) malloc(sizeof(ILXstatus)))==NULL)
    IndiLineExit("IndiLineExtractInit: memory allocation error");

  ILX->debugLevel = debugLevel;
  ILX->frameLen = frameLen;
  ILX->fSample = fSample;
  ILX->maxLineFreq = min(maxLineFreq,fSample*FMAX_FSAMPLE_RATIO);;
  ILX->maxLineAmpl = maxLineAmpl;
  ILX->maxNumLine = maxNumLine;
  ILX->maxNumEnv = maxNumEnv;
  ILX->maxNumHarm = maxNumHarm;
  ILX->maxNumHarmLine = maxNumHarmLine;
  ILX->maxNumNoisePara = maxNumNoisePara;
  ILX->maxNoiseFreq = maxNoiseFreq;
  ILX->maxNumAllLine = maxNumLine+maxNumHarm*maxNumHarmLine;

  if (ILX->debugLevel)
    printf("frmlen=%d fs=%9.3f maxf=%9.3f maxa=%9.3f numlin=%d\n",
	   frameLen,fSample,ILX->maxLineFreq,maxLineAmpl,maxNumLine);

  ILX->winLen = winLen = 2*frameLen;
  ILX->fftLen = 1;				/* transform length  */
  while (ILX->fftLen<winLen)
    ILX->fftLen*=2;
						/*  win. for const freq est. */

  ILX->addStart = (winLen-frameLen)/2;
  ILX->bufLen = ILX->fftLen;		    /* buffer length */

  if (ILX->debugLevel)
    printf("buflen=%d,addStart=%d\n",ILX->bufLen,ILX->addStart);

  /* allocate memory for arrays */

  if (
      /* fixed content arrays ("static") */
      (ILX->costab=(float*)malloc(2*ILX->fftLen*sizeof(float)))==NULL ||
      (ILX->longWin=(float*)malloc(ILX->winLen*sizeof(float)))==NULL ||
      /* arrays for return values */
      (ILX->envPara=(float**)malloc(ILX->maxNumEnv*sizeof(float*)))==NULL ||
      (ILX->envPara[0]=(float*)malloc(ILX->maxNumEnv*ENVPARANUM*
				      sizeof(float)))==NULL ||
      (ILX->lineFreq=(float*)malloc(ILX->maxNumAllLine*sizeof(float)))==NULL ||
      (ILX->lineAmpl=(float*)malloc(ILX->maxNumAllLine*sizeof(float)))==NULL ||
      (ILX->linePhase=(float*)malloc(ILX->maxNumAllLine*
				     sizeof(float)))==NULL ||
      (ILX->lineEnv=(int*)malloc(ILX->maxNumAllLine*sizeof(int)))==NULL ||
      (ILX->linePred=(int*)malloc(ILX->maxNumAllLine*sizeof(int)))==NULL ||
      (ILX->numHarmLine=(int*)malloc(ILX->maxNumHarm*sizeof(int)))==NULL ||
      (ILX->harmFreq=(float*)malloc(ILX->maxNumHarm*sizeof(float)))==NULL ||
      (ILX->harmFreqStretch=(float*)malloc(ILX->maxNumHarm*
					   sizeof(float)))==NULL ||
      (ILX->harmLineIdx=(int*)malloc(ILX->maxNumHarm*sizeof(int)))==NULL ||
      (ILX->harmEnv=(int*)malloc(ILX->maxNumHarm*sizeof(int)))==NULL ||
      (ILX->harmPred=(int*)malloc(ILX->maxNumHarm*sizeof(int)))==NULL ||
      (ILX->harmRate=(float*)malloc(ILX->maxNumHarm*sizeof(float)))==NULL ||
      (ILX->noisePara=(float*)malloc(ILX->maxNumNoisePara*
				     sizeof(float)))==NULL ||
      /* variable content arrays */
      (ILX->cos=(float*)malloc(ILX->bufLen*sizeof(float)))==NULL ||
      (ILX->sin=(float*)malloc(ILX->bufLen*sizeof(float)))==NULL ||
      (ILX->envWin=(float*)malloc(ILX->winLen*sizeof(float)))==NULL ||
      (ILX->x=(float*)malloc(ILX->bufLen*sizeof(float)))==NULL ||
      (ILX->xOrgRe=(float*)malloc(ILX->fftLen*sizeof(float)))==NULL ||
      (ILX->xOrgIm=(float*)malloc(ILX->fftLen*sizeof(float)))==NULL ||
      (ILX->xmodRe=(float*)malloc(ILX->bufLen*sizeof(float)))==NULL ||
      (ILX->xmodIm=(float*)malloc(ILX->bufLen*sizeof(float)))==NULL ||
      (ILX->xsyn=(float*)malloc(ILX->bufLen*sizeof(float)))==NULL ||
      (ILX->xsynRe=(float*)malloc(ILX->fftLen*sizeof(float)))==NULL ||
      (ILX->xsynIm=(float*)malloc(ILX->fftLen*sizeof(float)))==NULL ||
      (ILX->lineDf=(float*)malloc(ILX->maxNumAllLine*sizeof(float)))==NULL ||
      (ILX->lflag=(int *)malloc(ILX->fftLen*sizeof(int)))==NULL ||
      (ILX->prevLineCont=(int *)malloc(ILX->maxNumAllLine*
				       sizeof(int)))==NULL ||
      (ILX->harmCosFact=(float*)malloc(ILX->maxNumAllLine*
				       sizeof(float)))==NULL ||
      (ILX->harmSinFact=(float*)malloc(ILX->maxNumAllLine*
				       sizeof(float)))==NULL ||
      (ILX->harmFc=(float*)malloc(ILX->maxNumAllLine*sizeof(float)))==NULL ||
      (ILX->harmFe=(float*)malloc(ILX->maxNumAllLine*sizeof(float)))==NULL ||
      (ILX->harmPsy=(int*)malloc(ILX->maxNumAllLine*sizeof(int)))==NULL ||
      /* arrays for frame-to-frame memory */
      (ILX->prevEnvPara=(float**)malloc(ILX->maxNumEnv*
					sizeof(float*)))==NULL ||
      (ILX->prevEnvPara[0]=(float*)malloc(ILX->maxNumEnv*ENVPARANUM*
					  sizeof(float)))==NULL ||
      (ILX->prevLineFreq=(float*)malloc(ILX->maxNumAllLine*
					sizeof(float)))==NULL ||
      (ILX->prevLineAmpl=(float*)malloc(ILX->maxNumAllLine*
					sizeof(float)))==NULL ||
      (ILX->prevLinePhase=(float*)malloc(ILX->maxNumAllLine*
					 sizeof(float)))==NULL ||
      (ILX->prevLineEnv=(int*)malloc(ILX->maxNumAllLine*sizeof(int)))==NULL ||
      (ILX->prevNumHarmLine=(int*)malloc(ILX->maxNumHarm*sizeof(int)))==NULL ||
      (ILX->prevHarmFreq=(float*)malloc(ILX->maxNumHarm*
					sizeof(float)))==NULL ||
      (ILX->prevHarmFreqStretch=(float*)malloc(ILX->maxNumHarm*
					   sizeof(float)))==NULL ||
      (ILX->prevHarmLineIdx=(int*)malloc(ILX->maxNumHarm*sizeof(int)))==NULL ||
      (ILX->prevHarmEnv=(int*)malloc(ILX->maxNumHarm*sizeof(int)))==NULL)
    IndiLineExit("IndiLineExtractInit: memory allocation error");

  for (i=1;i<ILX->maxNumEnv;i++) {
    ILX->envPara[i] = ILX->envPara[0]+i*ENVPARANUM;
    ILX->prevEnvPara[i] = ILX->prevEnvPara[0]+i*ENVPARANUM;
  }
  
  /* init cosine-table */
  pf = M_PI/ILX->fftLen;
  for(i=0;i<2*ILX->fftLen;i++)
    ILX->costab[i] = cos(pf*i);

  /* init window functions and tables */
  for(i=0;i<winLen;i++)			/* long window */
    ILX->longWin[i] = .5*(1.-cos(M_PI*(2.*i+1.)/winLen));
  /* HP 980216   this could caused a seg fault ... thanks to Matti Koskinen
  for(i=winLen;i<ILX->fftLen;i++)
    ILX->longWin[i] = 0;
  */

  /* init frame-to-frame memory */
  ILX->prevNumLine = 0;

  ILX->prevHarmFreq[0] = 0;
  ILX->prevFundaF = 0;

  /* init envelope, line, harmonic tone parameter estimation */
  ILX->envStatus = ILXEnvInit(frameLen,&addStartEnv,&addEndEnv,
				   ILX->debugLevel);
  ILX->lineStatus = ILXLineInit(0);
  ILX->harmStatus = ILXHarmInit((ILX->debugLevel>>8)&3,
				ILX->fftLen,ILX->fSample);

  /* init psychoacoustic model */
  imax_psy = (int)(min(1.2*floor(ILX->fftLen*ILX->maxLineFreq/ILX->fSample),
			ILX->fftLen/2));

  UHD_init_psycho_acoustic(ILX->fSample,imax_psy*ILX->fSample/ILX->fftLen,
			   ILX->fftLen,frameLen,imax_psy,frameLen);

  return(ILX);
}


/* IndiLineExtract */
/* Extraction of individual lines and envelope from the current */
/* residual signal. */

void IndiLineExtract (
  ILXstatus *ILX,		/* in: ILX status handle */
  float *residualSignal,	/* in: residual signal */
				/*  [0..2.5*frameLen-1] */
				/*  previous frame, */
				/*  current frame, */
				/*  1st half next frame */
  float *synthSignal,		/* in: synthesised signal */
				/*  [0..2.5*frameLen-1] */
				/*  (currently not used) */
  float *finePitchPeriod,	/* in: (currently not used) */
  float *slowHarm,		/* in: (currently not used) */
  float *LPCcoef,		/* in: (currently not used) */
  int numLine,			/* in: num lines to extract */
				/*     -1 to disable extraction */
  int *numEnv,			/* out: num envelopes */
  float ***envPara,		/* out: envelope parameters */
				/*      [0..numEnv-1][0..ENVPARANUM-1] */
  int *numLineExtract,		/* out: num lines extracted */
  float **lineFreq,		/* out: line frequency [Hz] */
				/*      [0..numLine-1] */
  float **lineAmpl,		/* out: line amplitude */
				/*      [0..numLine-1] */
  float **linePhase,		/* out: line phase [rad] */
				/*      [0..numLine-1] */
  int **lineEnv,		/* out: line envelope flag/idx */
				/*      [0..numLine-1] */
  int **linePred,		/* out: line predecessor idx */
				/*      [0..numLine-1] */
  int *numHarm,			/* out: num harmonic tones */
  int **numHarmLine,		/* out: num lines of harmonic tone */
				/*      [0..numHarm-1] */
  float **harmFreq,		/* out: harm fundamental frequency [Hz] */
				/*      [0..numHarm-1] */
  float **harmFreqStretch,	/* out: harm freq stretching */
				/*      [0..numHarm-1] */
  int **harmLineIdx,		/* out: harm line para start idx */
				/*      [0..numHarm-1] */
				/*      ampl: lineAmpl[idx] etc. */
				/*      idx=harmLineIdx+(0..numHarmLine-1) */
  int **harmEnv,		/* out: harm envelope flag/idx */
				/*      [0..numHarm-1] */
  int **harmPred,		/* out: harm tone predecessor idx */
				/*      [0..numHarm-1] */
  float **harmRate,		/* out: bitrate for harm data (0..1) */
				/*      [0..numHarm-1] */
  int *numNoisePara,		/* out: num noise parameter */
  float *noiseFreq,		/* out: max noise freq (bandwidth) [Hz] */
  float **noisePara,		/* out: noise parameter (DCT) */
				/*      [0..numNoisePara-1] */
  int *noiseEnv,		/* out: noise envelope flag/idx */
  float *noiseRate)		/* out: bitrate for noise data (0..1) */

{
  int i,i0,iline,imin,imax,envFlag;	/* counters and indices */

  int frameCenter;			/* index of frame center */
  int constOffsC;			/* offset center window */

  float qx;				/* square sum of residual */

  /* actual line parameters */
  float actLineFreq,actLineAmpl,actLinePhase,actLineDf;
  int actLineEnv;


  float *win;

  /* variables for psychoacoustic model */
  float err;		/* spectral residual error */
  float emr,emrMax;	/* psychoacoustic criterion */
  float f0;		/* freq. of line with max. emr */
  float df0,fc,fe;	/* sweep rate, center & end frequency */
  float phi1,phi2,phiT;	/* for synthesis */

  /* parameters for lines with const. frequency  */
  float df,f,phi;
  float cfCosFact,cfSinFact,cfQeDiff;
  float cfCosFactEnv,cfSinFactEnv,cfQeDiffEnv;
  float qe;
  int cfEnvFlag;

  /* parameters for synthesis  */
  float *optCos,*optSin,optCosFact,optSinFact,xsynt,asn,aen;

  /* variables for line continuation */
  int bestPred;
  float envAmplE,envAmplS,as,ae,fd,fdq,dfq,fp,fa,q1,q2,q,bestQ;

  /* fundamental frequency */
  float fundaF,fundaRelSweep;
  float fundaFstretch;
  float fundaF_t;
  int iharm;
  int totline;
  float freq;
  float qxOri,qeSum,qeSumEnv;
  int numHarmPsy;
  int harmEnable,checkHarm;
  int di,j;
  int baseIdx;
  int addToMask;

  int harmIndi[MAXHARMINDI];

  int lineGood;
  int numLineBad;
  float harmPwr,harmPwr1,harmPwr2,harmPwr3;
  int harmCorr;
  float minHarmAmpl;

  float a,b,det,sfi1,sfi2,si2,si3,si4;

  /* spectrum pwr */
  float winNorm,specPwr,linePwr,envPwr;

  /* noise */
  float lfNoiseMean;
  float noiseEnvPwrFact;
  int noiseBW;

  /* init internal variables */
  frameCenter = ILX->addStart+ILX->frameLen/2;	/* index of frame center */
  constOffsC = frameCenter-ILX->winLen/2;

  /* index for minimum and maximum frequency */
  imin = min((int)(ILX->fftLen*ILX_FMIN/ILX->fSample),ILX->fftLen/2);
  imax = min((int)(ILX->fftLen*ILX->maxLineFreq/ILX->fSample),ILX->fftLen/2);

  /* core */

  if (ILX->debugLevel)
    printf("IndiLineExtract core\n");

  /* copy input samples */
  for (i=0; i<ILX->bufLen; i++)
    ILX->x[i]=(ILX->frameLen-ILX->addStart+i<5*ILX->frameLen/2
	       && ILX->frameLen-ILX->addStart+i>0)
      ? residualSignal[ILX->frameLen-ILX->addStart+i] : 0;
  /* current frame: x[addStart .. addStart+frameLen-1] */

  if (numLine==-1) {
    /* don't extract anything - return pointers to alloc'ed arrays */

    ILX->prevNumLine = 0;
    ILX->prevNumHarm = 0;
    
    *envPara = ILX->envPara;
    *lineFreq = ILX->lineFreq;
    *lineAmpl = ILX->lineAmpl;
    *linePhase = ILX->linePhase;
    *lineEnv = ILX->lineEnv;
    *linePred = ILX->linePred;
    *numHarmLine = ILX->numHarmLine;
    *harmFreq = ILX->harmFreq;
    *harmFreqStretch = ILX->harmFreqStretch;
    *harmLineIdx = ILX->harmLineIdx;
    *harmEnv = ILX->harmEnv;
    *harmPred = ILX->harmPred;
    *harmRate = ILX->harmRate;
    *noisePara = ILX->noisePara;
    
    return;
  }

  if (ILX->maxNumEnv == 0) {
    *numEnv = 0;
    envFlag = 0;
  }
  else {
    *numEnv = 1;
    /* calculate envelope parameters */
    envFlag = ILXEnvEstim(ILX->envStatus,&ILX->x[ILX->addStart],
			  ILX->envPara[0],ILX->envWin,ILX->winLen);

    if (ILX->debugLevel)
      printf("env: tm=%7.5f atk=%7.5f dec=%7.5f\n",
	     ILX->envPara[0][0],ILX->envPara[0][1],ILX->envPara[0][2]);
  }

  win = ILX->longWin;

  /* reset synth. signal */  
  memset((char *)ILX->xsyn, 0, ILX->bufLen*sizeof(float));
  
  winNorm = 0;
  for(i=0;i<ILX->winLen;i++)
    winNorm += win[i]*win[i];
  /* hanning: winNorm = 3./8.*winLen */

  qx = 0;
  for(i=0;i<ILX->winLen;i++)			/* calc. variance */
    qx += win[i]*win[i]*ILX->x[constOffsC+i]*ILX->x[constOffsC+i];
  qxOri = qx;
  if (ILX->debugLevel)
    printf(" qx=%e\n",qx);

  /* transform original signal to frequency domain */
  WinOdft(&ILX->x[constOffsC],win,ILX->winLen,ILX->fftLen,ILX->xOrgRe,ILX->xOrgIm);

  /* extract harmonic components ??? */
  if (ILX->maxNumHarm == 0) {
    *numHarm = 0;
    ILX->numHarmLine[0] = 0;
  }
  else {
    /* fundamental frequency estimation */
    fundaF = ILXHarmFundFreq(ILX->harmStatus,ILX->xOrgRe,ILX->xOrgIm);
    fundaFstretch = 0;
    fundaF_t = fundaF;
    
    fundaRelSweep = 0;
    if (fundaF && ILX->prevFundaF &&
	max(fundaF/ILX->prevFundaF,ILX->prevFundaF/fundaF) < 1.15)
      fundaRelSweep = fundaF/ILX->prevFundaF - 1.;
    ILX->prevFundaF = fundaF;
    if (ILX->debugLevel)
      printf("fundaF %7.1f fundaFstretch %8.5f fundaRelSweep %7.4f\n",
	     fundaF,fundaFstretch,fundaRelSweep);

    for (i=0; i<ILX->maxNumHarm; i++)
      ILX->harmLineIdx[i] = ILX->maxNumLine+i*ILX->maxNumHarmLine;

    *numHarm = 1;
    baseIdx = ILX->harmLineIdx[0];
    ILX->numHarmLine[0] = min(ILX->maxNumHarmLine,
			      (int)(min(ILX->fSample*FMAX_FSAMPLE_RATIO,
					ILX->maxLineFreq)/fundaF));
    if (ILX->debugLevel)
      printf("numHarmLine %2d\n",ILX->numHarmLine[0]);
  }

  iline = -1;
  qeSum = qeSumEnv = 0;
  harmPwr = harmPwr1 = harmPwr2 = harmPwr3 = 0.;

  /* harmonic line loop */
  for (iharm=0; iharm<ILX->numHarmLine[0]; iharm++) {

    /* transform synth. signal to frequency domain */
    WinOdft(&ILX->xsyn[constOffsC],win,ILX->winLen,ILX->fftLen,
	    ILX->xsynRe,ILX->xsynIm);

    f0 = (iharm+1)*(1+(iharm+1)*(iharm+1)*fundaFstretch)*
      2.*M_PI*(fundaF/ILX->fSample);
    df0 = f0*fundaRelSweep;
    i0 = (int)(f0/2./M_PI*ILX->fftLen);
    if (ILX->debugLevel)
      printf("harm idx %2d  harm freq %7.1f\n",
	     iharm+1,(iharm+1)*(1+(iharm+1)*(iharm+1)*fundaFstretch)*fundaF);

    specPwr = 0;
    for (i=max(0,i0-2); i<=min(ILX->fftLen/2-1,i0+2); i++)
      specPwr += ILX->xOrgRe[i]*ILX->xOrgRe[i]+ILX->xOrgIm[i]*ILX->xOrgIm[i];
    specPwr *= (2/(ILX->fftLen*winNorm));

    fc = f0;						/* center frequency */
    fe = f0+df0/2;					/* end frequency */


    /* limit fc & fe to FMIN..FMAX !!!   HP 960502 */
    fc = max(min(fc,ILX->maxLineFreq/ILX->fSample*2*M_PI),
		ILX_FMIN/ILX->fSample*2*M_PI);
    fe = max(min(fe,ILX->maxLineFreq/ILX->fSample*2*M_PI),
		ILX_FMIN/ILX->fSample*2*M_PI);

    /* calc. amplitudes and phases  - frequency sweep */
    phi1 = fc;
    /* phi1+2*phi2*ILX->frameLen/2 = fe */
    phi2 = (fe-phi1)/ILX->frameLen;
    /*      printf("fs=%7.2f, ",dphi*ILX->fSample/2/M_PI); */
    for(i=0;i<ILX->bufLen;i++) {		/* gen. complex sine wave */
      phiT = i-frameCenter+0.5;			/*  with frequency sweep */
      phi = phi1*phiT+phi2*phiT*phiT;
      ILX->cos[i] = cos(phi);
      ILX->sin[i] = sin(phi);
    }
    cfEnvFlag = envFlag && fc>ENV_MINFREQ*2*M_PI/ILX->fSample;
    CalcFactors(&ILX->x[constOffsC],&ILX->cos[constOffsC],
		&ILX->sin[constOffsC],win,
		cfEnvFlag ? ILX->envWin : (float*)NULL,
		ILX->winLen,
		&cfCosFact,&cfSinFact,&cfQeDiff,
		&cfCosFactEnv,&cfSinFactEnv,&cfQeDiffEnv,
		ILX->debugLevel);
    qeSum -= cfQeDiff;
    qeSumEnv -= cfEnvFlag ? cfQeDiffEnv : cfQeDiff;
    if (cfEnvFlag && cfQeDiffEnv < cfQeDiff) {
      cfCosFact = cfCosFactEnv;
      cfSinFact = cfSinFactEnv;
      cfQeDiff = cfQeDiffEnv;
    }
    else
      cfEnvFlag = 0;
    qe = qx + cfQeDiff;


    /* store parameters for best synthesis method */
    optCos = ILX->cos;
    optSin = ILX->sin;
    optCosFact = cfCosFact;
    optSinFact = cfSinFact;
    actLineFreq = fc*ILX->fSample/2/M_PI;		/* frequency in Hz */
    actLineDf = (fe-fc)*ILX->fSample/M_PI;		/* frequency change */
    actLineEnv = cfEnvFlag;
    qx = qe;					/* square sum residual error */
    actLineAmpl = sqrt(optCosFact*optCosFact+optSinFact*optSinFact);
    actLinePhase = atan2m(optSinFact,optCosFact);

    ILX->harmCosFact[baseIdx+iharm] = optCosFact;
    ILX->harmSinFact[baseIdx+iharm] = optSinFact;
    ILX->harmFc[baseIdx+iharm] = fc;
    ILX->harmFe[baseIdx+iharm] = fe;
    ILX->lineEnv[baseIdx+iharm] = cfEnvFlag;

    linePwr = actLineAmpl*actLineAmpl/2;
    if (actLineEnv)
      envPwr = EnvPwrFact(ILX->envPara[0]);
    else
      envPwr = 1;
    linePwr *= envPwr;
    if (ILX->debugLevel)
      printf("spec=%f line=%f l/s=%f env=%f\n",
	     specPwr,linePwr,linePwr/specPwr,envPwr);

    harmPwr += linePwr;
    if (linePwr/specPwr > HARMSPECTHRES) {
      harmPwr1 += linePwr;
      if ((iharm+1)%2 == 0)
	harmPwr2 += linePwr;
      if ((iharm+1)%3 == 0)
	harmPwr3 += linePwr;
    }

    /* remove spectral line */
    if(!actLineEnv) {			/* line with const. ampl. */
      for(i=0;i<ILX->bufLen;i++) {
	xsynt = optCosFact*optCos[i]-optSinFact*optSin[i];
	ILX->x[i] -= xsynt;
	ILX->xsyn[i] += xsynt;
      }
    }
    else {				/* line with ampl. envelope */
      asn = ILX->envWin[0];		/* start factor of env window */
      for(i=0;i<constOffsC;i++) {
	xsynt = asn*(optCosFact*optCos[i]-optSinFact*optSin[i]);
	ILX->x[i] -= xsynt;
	ILX->xsyn[i] += xsynt;
      }
      for(i=constOffsC;i<constOffsC+ILX->winLen;i++) {
	xsynt = ILX->envWin[-constOffsC+i]
	  *(optCosFact*optCos[i]-optSinFact*optSin[i]);
	ILX->x[i] -= xsynt;
	ILX->xsyn[i] += xsynt;
      }
      aen = ILX->envWin[ILX->winLen-1];	/* end factor of env window */
      for(i=constOffsC+ILX->winLen;i<ILX->bufLen;i++) {
	xsynt = aen*(optCosFact*optCos[i]-optSinFact*optSin[i]);
	ILX->x[i] -= xsynt;
	ILX->xsyn[i] += xsynt;
      }
    }

    if (ILX->debugLevel) {
      printf("h %2d ",iharm+1);
      printf("%s ",cfEnvFlag?" env ":"noenv");
      printf("sweep: f=%7.1f df=%7.1f a=%7.1f p=%5.2f e2=%7.1e\n",
	     fc*ILX->fSample/2/M_PI,
	     (fe-fc)*ILX->fSample/M_PI,
	     sqrt(cfCosFact*cfCosFact+cfSinFact*cfSinFact),
	     atan2m(cfSinFact,cfCosFact),qe);
    }

    /* store harmonic components */
    ILX->lineFreq[baseIdx+iharm] = actLineFreq;
    ILX->lineAmpl[baseIdx+iharm] = actLineAmpl;
      
  }
  /* end of harmonic line loop */
  
  /* improve funda freq estimate (*1, *2, *3) */
  harmCorr = 0;
  if (harmPwr1/harmPwr > HARMPWRTHRES) {
    harmCorr = 1;
    if (harmPwr2/harmPwr1 > HARMCORRTHRES)
      harmCorr = 2;
    if (harmPwr3/harmPwr1 > HARMCORRTHRES)
      harmCorr = 3;
  }

  if (ILX->debugLevel)
    printf("harmPwr (corr=%1d)  t=%10.2e  1=%10.2e  2=%10.2e  3=%10.2e\n",
	   harmCorr,harmPwr,harmPwr1,harmPwr2,harmPwr3);

  /* funda freq correction */
  if (harmCorr==0)
    ILX->numHarmLine[0] = 0;
  else
    ILX->numHarmLine[0] = ILX->numHarmLine[0]/harmCorr;
  if (harmCorr > 1)
    for (i=0; i<ILX->numHarmLine[0]; i++) {
      ILX->harmCosFact[baseIdx+i] = ILX->harmCosFact[baseIdx+(i+1)*harmCorr-1];
      ILX->harmSinFact[baseIdx+i] = ILX->harmSinFact[baseIdx+(i+1)*harmCorr-1];
      ILX->harmFc[baseIdx+i] = ILX->harmFc[baseIdx+(i+1)*harmCorr-1];
      ILX->harmFe[baseIdx+i] = ILX->harmFe[baseIdx+(i+1)*harmCorr-1];
      ILX->lineEnv[baseIdx+i] = ILX->lineEnv[baseIdx+(i+1)*harmCorr-1];
      ILX->lineFreq[baseIdx+i] = ILX->lineFreq[baseIdx+(i+1)*harmCorr-1];
      ILX->lineAmpl[baseIdx+i] = ILX->lineAmpl[baseIdx+(i+1)*harmCorr-1];
    }

  /* check for low-ampl harm lines */
  minHarmAmpl = 0;
  for (iharm=0; iharm<ILX->numHarmLine[0]; iharm++)
    if (ILX->lineAmpl[baseIdx+iharm] > minHarmAmpl)
      minHarmAmpl = ILX->lineAmpl[baseIdx+iharm];
  minHarmAmpl = max(ILX->maxLineAmpl*ILX_MINAMPL,minHarmAmpl*ILX_HARMMINAMPL);
  numHarmPsy = 0;
  for (iharm=0; iharm<ILX->numHarmLine[0]; iharm++)
    numHarmPsy += (ILX->lineAmpl[baseIdx+iharm] >= minHarmAmpl) ? 1 : 0;
  if (ILX->debugLevel)
    printf("numHarm gt minAmpl %2d\n",numHarmPsy);
  if (numHarmPsy<ILX_MINNUMHARM)
    ILX->numHarmLine[0] = 0;

  ILX->harmFreq[0] = 0;
  ILX->harmFreqStretch[0] = 0;
  if (ILX->numHarmLine[0]) {
    /* stretching: f(i) = f*i(1+s*i) = f*i + s*f*i^2   a=f   b=s*f */
    si2 = si3 = si4 = sfi1 = sfi2 = 0;
    for (iharm=0; iharm<ILX->numHarmLine[0]; iharm++)
      if (ILX->lineAmpl[baseIdx+iharm] >= minHarmAmpl) {
	/* HP 970902   f0,stretch only froum "loud" lines */
	i = iharm+1;
	si2 += (float)i*i;
	si3 += (float)i*i*i;
	si4 += (float)i*i*i*i;
	sfi1 += (float)i*ILX->lineFreq[baseIdx+iharm];
	sfi2 += (float)i*i*ILX->lineFreq[baseIdx+iharm];
      }
    if (ILX->numHarmLine[0] > 1) {
      det = si2*si4-si3*si3;
      a = (sfi1*si4-sfi2*si3)/det;
      b = (si2*sfi2-si3*sfi1)/det;
      ILX->harmFreq[0] = a;
      ILX->harmFreqStretch[0] = b/a;
    }
    else
      ILX->harmFreq[0] = ILX->lineFreq[baseIdx+0];
  }
  if (ILX->debugLevel)
    printf("harmFreq %7.1f  harmFreqStretch %8.1e\n",
	   ILX->harmFreq[0],ILX->harmFreqStretch[0]);
  if (ILX->debugLevel)
    for (i=1;i<=ILX->numHarmLine[0];i++)
      printf("stretch %2d %7.1f (f=%7.1f a=%7.1f)\n",
	     i-1,(ILX->harmFreq[0])*i*(1+(ILX->harmFreqStretch[0])*i),
	     ILX->lineFreq[baseIdx+i-1],ILX->lineAmpl[baseIdx+i-1]);

  ILX->harmEnv[0] = 0;	/* no envelope */
  if (ILX->numHarmLine[0] && envFlag && qeSumEnv>qeSum)
    ILX->harmEnv[0] = 1;
  ILX->harmPred[0] = 0;	/* new harm */
  if (ILX->prevHarmFreq[0]>0 && ILX->harmFreq[0]>0) {
    df = ILX->harmFreq[0]/ILX->prevHarmFreq[0];
    if (max(df,1/df)<1.15)
      ILX->harmPred[0] = 1;	/* continued harm */
  }
  if (ILX->debugLevel)
    printf("harmFreq %7.1f  harmEnv %d  harmPred %d  qe %7.1e  qeEnv %7.1e\n",
	   ILX->harmFreq[0],ILX->harmEnv[0],ILX->harmPred[0],qeSum,qeSumEnv);

  ILX->harmRate[0] = 0.;
  for (iharm=0; iharm<ILX->numHarmLine[0]; iharm++)
    ILX->harmPsy[baseIdx+iharm] = 0;

  /* copy input samples */
  for (i=0; i<ILX->bufLen; i++)
    ILX->x[i]=(ILX->frameLen-ILX->addStart+i<5*ILX->frameLen/2
	       && ILX->frameLen-ILX->addStart+i>0)
      ? residualSignal[ILX->frameLen-ILX->addStart+i] : 0;
  /* current frame: x[addStart .. addStart+frameLen-1] */

  /* reset synth. signal */  
  memset((char *)ILX->xsyn, 0, ILX->bufLen*sizeof(float));
  
  qx = qxOri;

  *numLineExtract = numLine;

  totline = 0;
  numHarmPsy = 0;
  numLineBad = 0;
  harmEnable = 0;
  addToMask = 0;
  /* individual line loop */
  for (iline=0; iline<numLine; iline++) {

    /* residual before first iteration = original */
    if (totline==0) {
      /* transform original signal to frequency domain */
      WinOdft(&ILX->x[constOffsC],win,ILX->winLen,ILX->fftLen,ILX->xOrgRe,ILX->xOrgIm);
    
      for (i=imin;i<=imax;i++)
	ILX->lflag[i] = 0;
    }

    /* transform synth. signal to frequency domain */
    WinOdft(&ILX->xsyn[constOffsC],win,ILX->winLen,ILX->fftLen,
	    ILX->xsynRe,ILX->xsynIm);

    /* calc. masking threshold caused by synth. signal (Re,Im,Chan) */
    ILX->excitSynthFft = UHD_psycho_acoustic(ILX->xsynRe,ILX->xsynIm,0); 

    emrMax = IL_MINFLOAT;
    i0 = 0;

    /* find line with maximum "error-to-mask-ratio" */
    for(i=imin;i<=imax;i++) {
      err = sqrt(ILX->xOrgRe[i]*ILX->xOrgRe[i]+ILX->xOrgIm[i]*ILX->xOrgIm[i])
	-sqrt(ILX->xsynRe[i]*ILX->xsynRe[i]+ILX->xsynIm[i]*ILX->xsynIm[i]);
      err = err*err;
      emr = err / ILX->excitSynthFft[i];	/* "error-to-mask-ratio" */

      if(!ILX->lflag[i] && emr>emrMax) {	/* max. value */
	i0 = i;
	emrMax = emr;
      }
    }

    freq = ILX->fSample*(i0+.5)/ILX->fftLen;
    iharm = -1;
    if (ILX->numHarmLine[0] > 0) {
      iharm = (int)(freq/(ILX->harmFreq[0])+.5)-1;
      iharm = (int)(freq/((ILX->harmFreq[0])*(1+(ILX->harmFreqStretch[0])*(iharm+1)))+.5)-1;
    }
    if (iharm >= ILX->numHarmLine[0])
      iharm = -1;

    if (ILX->debugLevel)
      printf("%2d: freq=%7.1f  ",totline,freq);

    specPwr = 0;
    for (i=max(0,i0-2); i<=min(ILX->fftLen/2-1,i0+2); i++)
      specPwr += ILX->xOrgRe[i]*ILX->xOrgRe[i]+ILX->xOrgIm[i]*ILX->xOrgIm[i];
    specPwr *= (2/(ILX->fftLen*winNorm));

    addToMask = 1;

    checkHarm = (iharm>=0 && ILX->harmPsy[baseIdx+iharm]==0 &&
		 fabs(freq-(ILX->harmFreq[0])*
		      (1+(ILX->harmFreqStretch[0])*(iharm+1))*(iharm+1))<20.);
    if (checkHarm) {
      /* take harmonic line */
      if (ILX->debugLevel)
	printf("harm %2d\n",iharm);

      ILX->harmPsy[baseIdx+iharm] = totline+1;

      if (numHarmPsy<MAXHARMINDI)
	harmIndi[numHarmPsy] = iline;
      if (numHarmPsy+1==MAXHARMINDI) {
	harmEnable = 1;
	/* remove MAXHARMINDI harmonic lines from indi line list */
	for (i=0; i<iline-MAXHARMINDI+1; i++) {
	  di = 0;
	  for (j=0; j<MAXHARMINDI; j++)
	    if (harmIndi[j] <= i+di)
	      di++;
	  if (di>0) {
	    ILX->lineFreq[i] = ILX->lineFreq[i+di];
	    ILX->lineAmpl[i] = ILX->lineAmpl[i+di];
	    ILX->linePhase[i] = ILX->linePhase[i+di];
	    ILX->lineEnv[i] = ILX->lineEnv[i+di];
	    ILX->lineDf[i] = ILX->lineDf[i+di];
	    if (ILX->debugLevel)
	      printf("moveindi %2d -> %2d\n",i+di,i);
	  }
	}
	iline -= MAXHARMINDI-1;
	if (ILX->debugLevel)
	  printf("set harmEnable  new iline=%2d (old-%d)\n",
		 iline,MAXHARMINDI-1);
      }
      numHarmPsy++;

      if (harmEnable) {
	fc = ILX->harmFc[baseIdx+iharm];
	fe = ILX->harmFe[baseIdx+iharm];
	phi1 = fc;
	/* phi1+2*phi2*ILX->frameLen/2 = fe */
	phi2 = (fe-phi1)/ILX->frameLen;
	/*      printf("fs=%7.2f, ",dphi*ILX->fSample/2/M_PI); */
	for(i=0;i<ILX->bufLen;i++) {		/* gen. complex sine wave */
	  phiT = i-frameCenter+0.5;		/*  with frequency sweep */
	  phi = phi1*phiT+phi2*phiT*phiT;
	  ILX->cos[i] = cos(phi);
	  ILX->sin[i] = sin(phi);
	}
	optCos = ILX->cos;
	optSin = ILX->sin;
	optCosFact = ILX->harmCosFact[baseIdx+iharm];
	optSinFact = ILX->harmSinFact[baseIdx+iharm];
	actLineEnv = ILX->lineEnv[baseIdx+iharm];
	actLineFreq = ILX->harmFreq[0]*(iharm+1);	/* frequency in Hz stretch???*/
	qx = qx;				/* should be reduced by qe */
      }
    }

    if (!checkHarm || !harmEnable) {
      /* estimate indiline parameter */
      iharm = -1;
      if (ILX->debugLevel)
	printf("indi %2d\n",iline);

      ILX->lflag[i0] = 1;			/* prevent further selection */

      f = 2*M_PI*(i0+.5)/ILX->fftLen;	/* frequency of corr. basis function */
      for(i=0;i<ILX->bufLen;i++) {	
	ILX->xmodRe[i] = ILX->x[i]*ILX->costab[((2*i0+1)*i)%(2*ILX->fftLen)];
	ILX->xmodIm[i] = ILX->x[i]
	  *ILX->costab[((2*i0+1)*i+ILX->fftLen/2)%(2*ILX->fftLen)];
      }

      /* limit f to FMIN..FMAX !!!   HP 960502 */
      f = max(min(f,ILX->maxLineFreq/ILX->fSample*2*M_PI),
	      ILX_FMIN/ILX->fSample*2*M_PI);

      /* calc. amplitudes and phases   - const. frequency */

      for(i=0;i<ILX->bufLen;i++) {		/* gen. complex sine wave */
	ILX->cos[i] = cos(f*(i-frameCenter+0.5));	/*  e**jf */
	ILX->sin[i] = sin(f*(i-frameCenter+0.5));
      }
      cfEnvFlag = envFlag && f>ENV_MINFREQ*2*M_PI/ILX->fSample;
      CalcFactors(&ILX->x[constOffsC],&ILX->cos[constOffsC],
		  &ILX->sin[constOffsC],win,
		  cfEnvFlag ? ILX->envWin : (float*)NULL,
		  ILX->winLen,
		  &cfCosFact,&cfSinFact,&cfQeDiff,
		  &cfCosFactEnv,&cfSinFactEnv,&cfQeDiffEnv,
		  ILX->debugLevel);
      if (cfEnvFlag && cfQeDiffEnv < cfQeDiff) {
	cfCosFact = cfCosFactEnv;
	cfSinFact = cfSinFactEnv;
	cfQeDiff = cfQeDiffEnv;
      }
      else
	cfEnvFlag = 0;
      qe = qx + cfQeDiff;

      /* store parameters for best synthesis method */
      optCos = ILX->cos;
      optSin = ILX->sin;
      optCosFact = cfCosFact;
      optSinFact = cfSinFact;
      actLineFreq = f*ILX->fSample/2/M_PI;		/* frequency in Hz */
      actLineDf = 0;
      actLineEnv = cfEnvFlag;
      qx = qe;					/* square sum residual error */

      if (ILX->debugLevel) {
	printf("%2d ",iline);
	printf("%s ",cfEnvFlag?" env ":"noenv");
	printf("const: f=%7.1f df=%7s a=%7.1f p=%5.2f e1=%7.1e\n",
	       f*ILX->fSample/2/M_PI,"",
	       sqrt(cfCosFact*cfCosFact+cfSinFact*cfSinFact),
	       atan2m(cfSinFact,cfCosFact),qe);
      }
    }

    actLineAmpl = sqrt(optCosFact*optCosFact+optSinFact*optSinFact);
    actLinePhase = atan2m(optSinFact,optCosFact);

    linePwr = actLineAmpl*actLineAmpl/2;
    if (actLineEnv)
      envPwr = EnvPwrFact(ILX->envPara[0]);
    else
      envPwr = 1;
    linePwr *= envPwr;
    lineGood = linePwr/specPwr > LINESPECTHRES;
    if (ILX->debugLevel)
      printf("lineGood %1d spec=%f line=%f l/s=%f env=%f\n",
	     lineGood,specPwr,linePwr,linePwr/specPwr,envPwr);

    /* remove spectral line */
    if(!actLineEnv) {			/* line with const. ampl. */
      for(i=0;i<ILX->bufLen;i++) {
	xsynt = optCosFact*optCos[i]-optSinFact*optSin[i];
	if (lineGood)
	  ILX->x[i] -= xsynt;
	if (addToMask)
	  ILX->xsyn[i] += xsynt;
      }
    }
    else {				/* line with ampl. envelope */
      asn = ILX->envWin[0];		/* start factor of env window */
      for(i=0;i<constOffsC;i++) {
	xsynt = asn*(optCosFact*optCos[i]-optSinFact*optSin[i]);
	if (lineGood)
	  ILX->x[i] -= xsynt;
	if (addToMask)
	  ILX->xsyn[i] += xsynt;
      }
      for(i=constOffsC;i<constOffsC+ILX->winLen;i++) {
	xsynt = ILX->envWin[-constOffsC+i]
	  *(optCosFact*optCos[i]-optSinFact*optSin[i]);
	if (lineGood)
	  ILX->x[i] -= xsynt;
	if (addToMask)
	  ILX->xsyn[i] += xsynt;
      }
      aen = ILX->envWin[ILX->winLen-1];	/* end factor of env window */
      for(i=constOffsC+ILX->winLen;i<ILX->bufLen;i++) {
	xsynt = aen*(optCosFact*optCos[i]-optSinFact*optSin[i]);
	if (lineGood)
	  ILX->x[i] -= xsynt;
	if (addToMask)
	  ILX->xsyn[i] += xsynt;
      }
    }

    totline++;
    if (iharm>=0)
      iline--;
    else {
      if (lineGood==0)
	numLineBad++;
      /* exit line loop if amplitude is too low */
      /* or if more than numline (!!!) "lineGood=0"-lines were encountered */
      if (actLineAmpl < ILX->maxLineAmpl*ILX_MINAMPL ||
	  numLineBad > numLine) {
	*numLineExtract = iline;
	if (ILX->debugLevel)
	  printf("exit line loop (ampl %8.2f < %8.2f)\n",
		 actLineAmpl,ILX->maxLineAmpl*ILX_MINAMPL);
	break;   /* exit line loop */
      }

      if (lineGood ||
	  (checkHarm && !harmEnable)) {	/* HP 970901, moveindi & lineGood=0 */
	/* output line parameters */
	ILX->lineFreq[iline] = actLineFreq;
	ILX->lineAmpl[iline] = actLineAmpl;
	ILX->linePhase[iline] = actLinePhase;
	ILX->lineEnv[iline] = actLineEnv;
	ILX->lineDf[iline] = actLineDf;
      }
      else
	iline--;
    }

  }
  /* end of individual line loop */

  /* check num harm used in il loop */
  ILX->harmRate[0] = numHarmPsy/(float)(numHarmPsy+*numLineExtract);
  if (ILX->debugLevel)
    printf("totline %2d  numLineBad %2d  numHarmPsy %2d  numLineExtract %2d\n"
	   "harmEnable %1d  harmRate %5.3f\n",
	   totline,numLineBad,numHarmPsy,*numLineExtract,
	   harmEnable,ILX->harmRate[0]);

  if (!harmEnable)
    ILX->numHarmLine[0] = 0;

  if (ILX->numHarmLine[0] == 0)
    *numHarm = 0;

  /* remove remaining harm lines from residual signal */
  if (ILX->debugLevel) {
    q = 0;
    for(i=0;i<ILX->winLen;i++)			/* calc. variance */
      q += win[i]*win[i]*ILX->x[constOffsC+i]*ILX->x[constOffsC+i];
    printf("remain harm   pre qx=%e",q);
  }

  if (harmEnable) {
    for (iharm=0; iharm<ILX->numHarmLine[0]; iharm++) {
      if (ILX->harmPsy[baseIdx+iharm]==0) {

	fc = ILX->harmFc[baseIdx+iharm];
	fe = ILX->harmFe[baseIdx+iharm];
	phi1 = fc;
	/* phi1+2*phi2*ILX->frameLen/2 = fe */
	phi2 = (fe-phi1)/ILX->frameLen;
	/*      printf("fs=%7.2f, ",dphi*ILX->fSample/2/M_PI); */
	for(i=0;i<ILX->bufLen;i++) {		/* gen. complex sine wave */
	  phiT = i-frameCenter+0.5;		/*  with frequency sweep */
	  phi = phi1*phiT+phi2*phiT*phiT;
	  ILX->cos[i] = cos(phi);
	  ILX->sin[i] = sin(phi);
	}

	actLineEnv = ILX->lineEnv[baseIdx+iharm];
	CalcFactors(&ILX->x[constOffsC],&ILX->cos[constOffsC],
		    &ILX->sin[constOffsC],win,
		    actLineEnv ? ILX->envWin : (float*)NULL,
		    ILX->winLen,
		    &cfCosFact,&cfSinFact,&cfQeDiff,
		    &cfCosFactEnv,&cfSinFactEnv,&cfQeDiffEnv,
		    ILX->debugLevel);

	optCos = ILX->cos;
	optSin = ILX->sin;

	if(!actLineEnv) {			/* line with const. ampl. */
	  optCosFact = cfCosFact;
	  optSinFact = cfSinFact;
	}
	else {
	  optCosFact = cfCosFactEnv;
	  optSinFact = cfSinFactEnv;
	}

	/* remove spectral line */
	if(!actLineEnv) {			/* line with const. ampl. */
	  for(i=0;i<ILX->bufLen;i++) {
	    xsynt = optCosFact*optCos[i]-optSinFact*optSin[i];
	    ILX->x[i] -= xsynt;
	  }
	}
	else {				/* line with ampl. envelope */
	  asn = ILX->envWin[0];		/* start factor of env window */
	  for(i=0;i<constOffsC;i++) {
	    xsynt = asn*(optCosFact*optCos[i]-optSinFact*optSin[i]);
	    ILX->x[i] -= xsynt;
	  }
	  for(i=constOffsC;i<constOffsC+ILX->winLen;i++) {
	    xsynt = ILX->envWin[-constOffsC+i]
	      *(optCosFact*optCos[i]-optSinFact*optSin[i]);
	    ILX->x[i] -= xsynt;
	  }
	  aen = ILX->envWin[ILX->winLen-1];	/* end factor of env window */
	  for(i=constOffsC+ILX->winLen;i<ILX->bufLen;i++) {
	    xsynt = aen*(optCosFact*optCos[i]-optSinFact*optSin[i]);
	    ILX->x[i] -= xsynt;
	  }
	}

      }
    }
  }

  if (ILX->debugLevel) {
    q = 0;
    for(i=0;i<ILX->winLen;i++)			/* calc. variance */
      q += win[i]*win[i]*ILX->x[constOffsC+i]*ILX->x[constOffsC+i];
    printf("   post qx=%e\n",q);
  }


#ifdef ASDFASDFASDF
  /* gen 1..5 harm at f0 ... */
  if (numLine <= 5) {
    fundaA = 0;
    for (iline=0; iline<numLine; iline++)
      if (ILX->lineAmpl[iline] > fundaA)
	fundaA = ILX->lineAmpl[iline];
    *numLineExtract = numLine;
    for (iline=0; iline<numLine; iline++) {
      ILX->lineFreq[iline] = (iline+1)*fundaF;
      ILX->lineAmpl[iline] = fundaA/(iline+1);
      ILX->lineAmpl[iline] = fundaA;
      ILX->linePhase[iline] = 0;
      ILX->lineEnv[iline] = 0;
      ILX->lineDf[iline] = 0;
    }
  }

  for (iline=0; iline<*numLineExtract; iline++)
    ILX->lineFreq[iline] = (iline+1)*fundaF_t;
#endif

  
  /* find predecessors for continued lines */
  
  envAmplE = max(MINENVAMPL,1-((1-ILX->prevEnvPara[0][0])*ILX->prevEnvPara[0][2]));
  envAmplS = max(MINENVAMPL,1-(ILX->envPara[0][0]*ILX->envPara[0][1]));
  for (i=0;i<ILX->prevNumLine;i++)
    ILX->prevLineCont[i]=0;
  for (iline=0;iline<*numLineExtract;iline++) {
    as = ILX->lineAmpl[iline]*(ILX->lineEnv[iline]?envAmplS:1);
    fdq = ILX->lineFreq[iline];
    fd = fdq-ILX->lineDf[iline];
    bestPred = 0;
    bestQ = 0;
    for (i=0;i<ILX->prevNumLine;i++) {
      if (!ILX->prevLineCont[i]) {
	ae = ILX->prevLineAmpl[i]*(ILX->prevLineEnv[i]?envAmplE:1);
	fp = ILX->prevLineFreq[i];
	/* calc q */
	if (ae>0 && as>0 && fp>0 && fd>0) {
	  dfq = fabs(fdq-fp)/fp;
	  df = fabs(fd-fp)/fp;
	  fa = as/ae;
	  if(dfq<DFQMAX && df<DFMAX && fa>FAMIN && fa<FAMAX) {
	    q1 = 1.-df/DFMAX;
	    if(fa<1)
	      q2 = (fa-FAMIN)/(1-FAMIN);
	    else
	      q2 = (1./fa-FAMIN)/(1-FAMIN);
	    q = q1*q2;
	    if (q > bestQ) {
	      bestQ = q;
	      bestPred = i+1;
	    }
	  }
	}
      }
    }
    ILX->linePred[iline] = bestPred;
    if (bestPred)
      ILX->prevLineCont[bestPred-1] = 1;
  }

  /* calc. parameters for residual noise */
  if (ILX->maxNumNoisePara == 0)
    *numNoisePara = 0;
  else {
    if (ILX->maxNumEnv <= 1) {
      *noiseEnv = 0;
      noiseEnvPwrFact = 1;
    }
    else {
      *numEnv = 2;
      /* calculate new envelope parameters for the residual signal */
      if (ILXEnvEstim(ILX->envStatus,&ILX->x[ILX->addStart],
		      ILX->envPara[1],ILX->envWin,ILX->winLen))
	*noiseEnv = 2;
      else
	*noiseEnv = 0;
#ifdef NOISE_ENV_THRESHOLD
      /* doesn't help, simply use me=1 (default) also for casta !!! */
      if (ILX->envPara[1][0]*ILX->envPara[1][1] < 1 &&
	  (1-ILX->envPara[1][0])*ILX->envPara[1][2] < 1) {
	/* env: a(t)>0 for all t=0..1 -> don't use this env */
	*noiseEnv = 0;
	ILX->envPara[1][0] = ILX->envPara[1][1] = ILX->envPara[1][2] = 0;
      }
#endif
      noiseEnvPwrFact = EnvPwrFact(ILX->envPara[1]);
    }
    
    /* residual variance */
    qx = 0;
    for(i=0;i<ILX->winLen;i++) {
      qx += win[i]*win[i]*ILX->x[constOffsC+i]*ILX->x[constOffsC+i];
    }

    WinOdft(&ILX->x[constOffsC],win,ILX->winLen,ILX->fftLen,
	    ILX->xsynRe,ILX->xsynIm);
    for(i=0;i<ILX->fftLen/2;i++)
      ILX->xsynRe[i] = sqrt(ILX->xsynRe[i]*ILX->xsynRe[i]+
			    +ILX->xsynIm[i]*ILX->xsynIm[i])/noiseEnvPwrFact;
    lfNoiseMean = 0;
    for(i=3;i<10;i++)
      lfNoiseMean += ILX->xsynRe[i];
    lfNoiseMean /= 8;
    ILX->xsynRe[0] = lfNoiseMean;
    ILX->xsynRe[1] = .7*lfNoiseMean+.3*ILX->xsynRe[1];
    ILX->xsynRe[2] = .3*lfNoiseMean+.7*ILX->xsynRe[2];
    
    noiseBW = (int)(ILX->fftLen/ILX->fSample*ILX->maxNoiseFreq+.5);
    NoiseDct(ILX->xsynRe,ILX->noisePara,min(ILX->fftLen/2,noiseBW),
	     ILX->maxNumNoisePara,noiseBW);
    
    for(i=0;i<ILX->maxNumNoisePara;i++)
      ILX->noisePara[i] *= 512.0/ILX->fftLen;
      /* HP 980216   compensate /=fftLen in ODFT (256 sample/frame)*/
    *numNoisePara = ILX->maxNumNoisePara;
    *noiseRate = qx/max(qxOri,1.);
    *noiseFreq = ILX->maxNoiseFreq;
    
    if (ILX->debugLevel) {
      printf("qxOri=%e, qx=%e, noiseRate=%4.2f\n",qxOri,qx,*noiseRate);
      /*     printf("NoisePara="); */
      /*     for(i=0;i<ILX->maxNoise;i++) */
      /*       printf(" %6.1f",ILX->noisePara[i]); */
      /*     printf("\n"); */
    }
  }

  /* set pointers to arrays with return values */

  *envPara = ILX->envPara;
  *lineFreq = ILX->lineFreq;
  *lineAmpl = ILX->lineAmpl;
  *linePhase = ILX->linePhase;
  *lineEnv = ILX->lineEnv;
  *linePred = ILX->linePred;
  *numHarmLine = ILX->numHarmLine;
  *harmFreq = ILX->harmFreq;
  *harmFreqStretch = ILX->harmFreqStretch;
  *harmLineIdx = ILX->harmLineIdx;
  *harmEnv = ILX->harmEnv;
  *harmPred = ILX->harmPred;
  *harmRate = ILX->harmRate;
  *noisePara = ILX->noisePara;

  /* copy parameters into frame-to-frame memory */

  ILX->prevNumLine = *numLineExtract;
  for (j=0; j<*numEnv; j++)
    for (i=0; i<ENVPARANUM; i++)
      ILX->prevEnvPara[j][i] = ILX->envPara[j][i];
  for (iline=0; iline<*numLineExtract; iline++) {
    ILX->prevLineFreq[iline] = ILX->lineFreq[iline];
    ILX->prevLineAmpl[iline] = ILX->lineAmpl[iline];
    ILX->prevLinePhase[iline] = ILX->linePhase[iline];
    ILX->prevLineEnv[iline] = ILX->lineEnv[iline];
  }
  ILX->prevNumHarm = *numHarm;
  for (i=0; i<*numHarm; i++) {
    ILX->prevNumHarmLine[i] = ILX->numHarmLine[i];
    ILX->prevHarmLineIdx[i] = ILX->harmLineIdx[i];
    ILX->prevHarmFreq[i] = ILX->harmFreq[i];
    ILX->prevHarmFreqStretch[i] = ILX->harmFreqStretch[i];
    ILX->prevHarmEnv[i] = ILX->harmEnv[i];
    for (iline=ILX->harmLineIdx[i];
	 iline<ILX->harmLineIdx[i]+ILX->numHarmLine[i]; iline++) {
      ILX->prevLineFreq[iline] = ILX->lineFreq[iline];
      ILX->prevLineAmpl[iline] = ILX->lineAmpl[iline];
      ILX->prevLinePhase[iline] = ILX->linePhase[iline];
      ILX->prevLineEnv[iline] = ILX->lineEnv[iline];
    }
  }
  

  if (ILX->debugLevel)
    for (i=0;i<*numLineExtract;i++)
      printf("%2d: f=%7.1f (df=%7.1f) a=%7.1f p=%5.2f e=%1d c=%2d\n",
	     i,ILX->lineFreq[i],ILX->lineDf[i],ILX->lineAmpl[i],
	     ILX->linePhase[i],ILX->lineEnv[i],ILX->linePred[i]);

}
  

/* IndiLineExtractFree() */
/* Free memory allocated by IndiLineExtractInit(). */

void IndiLineExtractFree (
  ILXstatus *ILX)		/* in: ILX status handle */
{
  free(ILX->costab);
  free(ILX->longWin);
  free(ILX->envPara[0]);
  free(ILX->envPara);
  free(ILX->lineFreq);
  free(ILX->lineAmpl);
  free(ILX->linePhase);
  free(ILX->lineEnv);
  free(ILX->linePred);
  free(ILX->numHarmLine);
  free(ILX->harmFreq);
  free(ILX->harmFreqStretch);
  free(ILX->harmLineIdx);
  free(ILX->harmEnv);
  free(ILX->harmPred);
  free(ILX->harmRate);
  free(ILX->noisePara);
  free(ILX->cos);
  free(ILX->sin);
  free(ILX->envWin);
  free(ILX->x);
  free(ILX->xOrgRe);
  free(ILX->xOrgIm);
  free(ILX->xmodRe);
  free(ILX->xmodIm);
  free(ILX->xsyn);
  free(ILX->xsynRe);
  free(ILX->xsynIm);
  free(ILX->lineDf);
  free(ILX->lflag);
  free(ILX->prevLineCont);
  free(ILX->harmCosFact);
  free(ILX->harmSinFact);
  free(ILX->harmFc);
  free(ILX->harmFe);
  free(ILX->harmPsy);
  free(ILX->prevEnvPara[0]);
  free(ILX->prevEnvPara);
  free(ILX->prevLineFreq);
  free(ILX->prevLineAmpl);
  free(ILX->prevLinePhase);
  free(ILX->prevLineEnv);
  free(ILX->prevNumHarmLine);
  free(ILX->prevHarmFreq);
  free(ILX->prevHarmFreqStretch);
  free(ILX->prevHarmLineIdx);
  free(ILX->prevHarmEnv);

  ILXEnvFree(ILX->envStatus);
  ILXLineFree(ILX->lineStatus);
  ILXHarmFree(ILX->harmStatus);

  free(ILX);
}

/* end of indilinextr.c */
