/*****************************************************************************/
/* This software module was originally developed by                          */
/*   Naoki Iwakami (NTT)                                                     */
/* and edited by                                                             */
/*   Naoki Iwakami (NTT) on 1997-07-17,                                      */
/* 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)1996.                                                        */
/*****************************************************************************/



#include "ntt_conf.h"
#include "ntt_enc_para.h"
#include "ntt_encode.h"

void ntt_tf_perceptual_model(/* Input */
			 double lpc_spectrum[],      /* LPC spectrum*/
			 double bark_env[],          /* Bark-scale envelope*/
			 double gain[],              /* gain factor */
			 int       w_type,           /* block type */
			 double spectrum[],          /* spectrum */
			 double pitch_sequence[],    /* periodic peak components */
			 /* Output */
			 double perceptual_weight[]) /* perceptual weight */
{
    /*--- Variables ---*/
    int    ismp, i_sup, top, isf, subtop;
    double ratio;
    double ratio0;

    switch(w_type){
    case ONLY_LONG_WINDOW:
    case LONG_SHORT_WINDOW:
    case SHORT_LONG_WINDOW:
    case LONG_MEDIUM_WINDOW:
    case MEDIUM_LONG_WINDOW:
	for (i_sup=0; i_sup<ntt_N_SUP; i_sup++){
	    top = i_sup * ntt_N_FR;
	    for (ismp=0; ismp<ntt_N_FR; ismp++){
		perceptual_weight[ismp+top] = 1./lpc_spectrum[ismp+top];
	    }
	    ntt_prcptw(perceptual_weight+top, bark_env+top);
            for (ismp=0; ismp<ntt_N_FR; ismp++){
		perceptual_weight[ismp+top] *=
		    lpc_spectrum[ismp+top]/bark_env[ismp+top];
	    }

	    for (ismp=0; ismp<ntt_N_FR; ismp++){
		ratio0= (spectrum[ismp+top]/lpc_spectrum[ismp+top]);
		  ratio = ((ratio0-pitch_sequence[ismp+top])*
		           (ratio0-pitch_sequence[ismp+top])+10.001)
		           /((ratio0*ratio0)+10.001);
		if(ratio<1.0) perceptual_weight[ismp+top] *= ratio;
	     }
	}
	break;
    case ONLY_MEDIUM_WINDOW:
    case MEDIUM_SHORT_WINDOW:
    case SHORT_MEDIUM_WINDOW:
	for (i_sup=0; i_sup<ntt_N_SUP; i_sup++){
	    top = i_sup * ntt_N_FR;
	    for (isf=0; isf<ntt_N_MID; isf++){
		subtop = top + isf * ntt_N_FR_M;
		for (ismp=0; ismp<ntt_N_FR_M; ismp++){
		    perceptual_weight[ismp+subtop] =
			bark_env[ismp+subtop] / lpc_spectrum[ismp+subtop];
		}
	    }
	}
	ntt_prcptw_m(perceptual_weight, gain, perceptual_weight);
	for (i_sup=0; i_sup<ntt_N_SUP; i_sup++){
	    top = i_sup * ntt_N_FR;
	    for (isf=0; isf<ntt_N_MID; isf++){
		subtop = top + isf * ntt_N_FR_M;
		for (ismp=0; ismp<ntt_N_FR_M; ismp++){
		    perceptual_weight[ismp+subtop] *=
			lpc_spectrum[ismp+subtop] / bark_env[ismp+subtop];
		}
	    }
	}
	break;
    case ONLY_SHORT_WINDOW:
	for (i_sup=0; i_sup<ntt_N_SUP; i_sup++){
	    top = i_sup * ntt_N_FR;
	    for (isf=0; isf<ntt_N_SHRT; isf++){
		subtop = top + isf * ntt_N_FR_S;
		for (ismp=0; ismp<ntt_N_FR_S; ismp++){
		    perceptual_weight[ismp+subtop] =
			bark_env[ismp+subtop] / lpc_spectrum[ismp+subtop];
		}
	    }
	}
	ntt_prcptw_s(perceptual_weight, gain, perceptual_weight);
	for (i_sup=0; i_sup<ntt_N_SUP; i_sup++){
	    top = i_sup * ntt_N_FR;
	    for (isf=0; isf<ntt_N_SHRT; isf++){
		subtop = top + isf * ntt_N_FR_S;
		for (ismp=0; ismp<ntt_N_FR_S; ismp++){
		    perceptual_weight[ismp+subtop] *=
			lpc_spectrum[ismp+subtop] / bark_env[ismp+subtop];
		}
	    }
	}
	break;
    }
    
}


void ntt_prcptw(double pwt[],   /* Input/Output -- Perceptual weighting factor */
	    double pred[])   /* Input        -- Interframe prediction factor */

{
    /*--- Variables ---*/
    int	ismp;
    double gamma_w, gamma_w_mic;

    gamma_w     = ntt_GAMMA_W;
    gamma_w_mic = ntt_GAMMA_W_MIC;

    /*--- Main operation ---*/
    for ( ismp=0; ismp<ntt_N_FR; ismp++ ){
      pwt[ismp] =  pow(pred[ismp], gamma_w) * pow(pwt[ismp], gamma_w_mic);
    }
}

void ntt_prcptw_m(double wt[],	/* Input  -- LPC weighting factor */
	      double gain[],	/* Input  -- Gain parameter */
	      double pwt[])	/* Output -- Perceptual weighting factor */
{
    /*--- Variables ---*/
    int		i_sup, ismp, i_mid, iptr, top;
    double	gtmp;

    /*--- Main operation ---*/
    for(i_sup=0; i_sup<ntt_N_SUP; i_sup++){
	for (i_mid=0; i_mid<ntt_N_MID; i_mid++){
	    iptr = i_sup * ntt_N_MID + i_mid;

	    gtmp =  pow(gain[iptr]+1.e-5,ntt_GAMMA_W_M_T) / (gain[iptr]+1.e-5);

	    top = iptr * ntt_N_FR_M;
	    for ( ismp=0; ismp<ntt_N_FR_M; ismp++ ){
               /* Temporal perceptual weighting is added */
		pwt[ismp+top] = gtmp*pow(wt[ismp+top], ntt_GAMMA_W_M); 
	    }
	}
    }
}

void ntt_prcptw_s(double wt[],	/* Input  -- LPC weighting factor */
	      double gain[],	/* Input  -- Gain parameter */
	      double pwt[])	/* Output -- Perceptual weighting factor */

{
    /*--- Variables ---*/
    int		i_sup, ismp, i_shrt, iptr, top;
    double	gtmp, gamma_w_s_t, gamma_w_s;

    gamma_w_s_t = ntt_GAMMA_W_S_T;
    gamma_w_s = ntt_GAMMA_W_S;

    /*--- Main operation ---*/
    for(i_sup=0; i_sup<ntt_N_SUP; i_sup++){
	for (i_shrt=0; i_shrt<ntt_N_SHRT; i_shrt++){
	    iptr = i_sup * ntt_N_SHRT + i_shrt;

	    gtmp =  pow(gain[iptr]+1.e-5,gamma_w_s_t) / (gain[iptr]+1.e-5);

	    top = iptr * ntt_N_FR_S;
	    for ( ismp=0; ismp<ntt_N_FR_S; ismp++ ){
               /* Temporal perceptual weighting is added */
		pwt[ismp+top] = gtmp*pow(wt[ismp+top], gamma_w_s);
	    }
	}
    }
}
