stft_window.c

LastWave

/*..........................................................................*/
/*                                                                          */
/*      L a s t W a v e    P a c k a g e 'stft' 2.1                         */
/*                                                                          */
/*      Copyright (C) 1997-2002 R.Gribonval, E.Bacry and J.Abadia           */
/*      email  : [e-mail]                                    */
/*               [e-mail]                            */
/*                                                                          */
/*..........................................................................*/
/*                                                                          */
/*      This program is a 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 (in a file named COPYRIGHT);                */
/*      if not, write to the Free Software Foundation, Inc.,                */
/*      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA             */
/*                                                                          */
/*..........................................................................*/


#include "lastwave.h"

#include "stft.h"


/*#define DEBUGSHIFT*/

/************************************************************/
/*     The windows                                          */
/*       the following functions must                       */
/*       fill up a signal                                   */
/*       corresponding to f(i+shift)                        */
/*       (starting i = 0)                                   */ 
/*       The signal must be normalized                      */
/*       The return value is the normalization factor       */
/************************************************************/


/*************************/
/*    Order 0 spline     */
/*************************/

static LWFLOAT Spline0(SIGNAL window,LWFLOAT shift)
{
  int j;
  LWFLOAT i;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Spline0 : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Spline0 : bad shift %g",shift);
  
  /* The window */
  for(j = 0, i=shift; i < window->size; j++, i += 1.0) {
    window->Y[j] = 1.0;
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero */
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  factor = 1/sqrt(energy);
  for(j = 0; j < window->size; j++) {
    window->Y[j] *= factor;
  }
  
  return(factor);
}

/*************************/
/*    Order 1 spline     */
/*************************/

static LWFLOAT Spline1(SIGNAL window,LWFLOAT shift)
{
  int j;
  LWFLOAT i;
  LWFLOAT temp;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Spline1 : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Spline1 : bad shift %g",shift);
  
  temp = window->size/2.0;
  /* The window */
  for(j=0, i = shift; i < window->size/2; j++, i += 1.0) {
    window->Y[j] = i/temp;
    energy += window->Y[j]*window->Y[j];
  }
  for(; i < window->size; j++, i += 1.0) {
    window->Y[j] = (2-i/temp);
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero */
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  factor = 1/sqrt(energy);
  for(j=0; j < window->size; j++) {
    window->Y[j] *= factor;
  }
  
  return(factor);
}

/*************************/
/*    Order 2 spline     */
/*************************/

static LWFLOAT Spline2(SIGNAL window,LWFLOAT shift)
{
  int j;
  LWFLOAT i;
  LWFLOAT temp;
  LWFLOAT x;
  LWFLOAT energy = 0;
  LWFLOAT factor;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Spline2 : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Spline2 : bad shift %g",shift);
  
  temp = window->size/3.0;
  
  /* The window */
  for(i = shift, j=0; i<temp; i++,j++) {
    x = i/temp;
    window->Y[j] = .5*x*x;
    energy += window->Y[j]*window->Y[j];
  }
  for(;i<2*temp;i++,j++) {
    x = i/temp-1.5;
    window->Y[j] = (3.0/4.0 - x*x);
    energy += window->Y[j]*window->Y[j];
  }
  for( ; i < window->size && j < window->size; i++,j++) {
    x = 3-i/temp;
    window->Y[j] = .5*x*x;
    energy += window->Y[j]*window->Y[j];
  }
  /* Putting the first point to zero */
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  
  /* Normalizing */
  factor = 1/sqrt(energy);
  for(j=0;j<window->size;j++) {
    window->Y[j] *= factor;
  }
  
  return(factor);
}


/*************************/
/*    Order 3 spline     */
/*************************/

static LWFLOAT Spline3(SIGNAL window,LWFLOAT shift)
{
  LWFLOAT i;
  int j;
  LWFLOAT temp = window->size/4.0;
  LWFLOAT x;
  LWFLOAT energy = 0;
  LWFLOAT factor;
  
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Spline3 : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Spline3 : bad shift %g",shift);
  
  /* The window */
  for(i=shift, j=0; i<temp; i++,j++) {
    x = i/temp;
    window->Y[j] = x*x*x/6;
    energy += window->Y[j]*window->Y[j];
  }
  for(; i < 2*temp; i++,j++) {
    x = i/temp;
    window->Y[j] = -.5*x*x*x+2*x*x-2*x+2.0/3;
    energy += window->Y[j]*window->Y[j];
  }
  for(; i < 3*temp; i++,j++) {
    x = 4-i/temp;
    window->Y[j] = -.5*x*x*x+2*x*x-2*x+2.0/3;
    energy += window->Y[j]*window->Y[j];
  }
  for(; i < window->size; i++,j++) {
    x = 4-i/temp;
    window->Y[j] = x*x*x/6;
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero */
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  factor = 1/sqrt(energy);
  for(j=0; j < window->size; j++) {
    window->Y[j] *= factor;
  }
  
  return(factor);
}


/*************************/
/*   Gaussian            */
/*************************/

/*double theSigma2 = 0.014;   */
double theGaussianSigma2 = 0.02;  

static LWFLOAT Gauss(SIGNAL window,LWFLOAT shift)
{
  LWFLOAT i;
  int   j;
  LWFLOAT x;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Gauss : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Gauss : bad shift %g",shift);
  
  /* The window */
  for(i = shift, j=0; j < window->size; j++,i++) {
    x = (i-window->size/2)/window->size;
    window->Y[j] = exp(-x*x/(2*theGaussianSigma2));
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero */
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  factor = 1/sqrt(energy);
  for(j=0;j<window->size;j++) {
    window->Y[j] *= factor;
  }
  return(factor);
}


static LWFLOAT Hamming(SIGNAL window,LWFLOAT shift)
{
  LWFLOAT i;
  int   j;
  LWFLOAT x;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Hamming : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Hamming : bad shift %g",shift);
  
  /* The window */
  for(i = shift, j=0; j < window->size; j++,i++) {
    x = (i-window->size/2)/window->size;
    window->Y[j] = 0.54+0.46*cos(2*M_PI*x);
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero */
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  factor = 1/sqrt(energy);
  for(j=0;j<window->size;j++) {
    window->Y[j] *= factor;
  }
  return(factor);
}

static LWFLOAT Hanning(SIGNAL window,LWFLOAT shift)
{
  LWFLOAT i;
  int   j;
  LWFLOAT x;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Hanning : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Hanning : bad shift %g",shift);
  
  /* The window */
  for(i = shift, j=0; j < window->size; j++,i++) {
    x = (i-window->size/2)/window->size;
    window->Y[j] = cos(M_PI*x);
    window->Y[j] *= window->Y[j];
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero */
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  factor = 1/sqrt(energy);
  for(j=0;j<window->size;j++) {
    window->Y[j] *= factor;
  }
  return(factor);
}

static LWFLOAT Blackman(SIGNAL window,LWFLOAT shift)
{
  LWFLOAT i;
  int   j;
  LWFLOAT x;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Blackman : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Blackman : bad shift %g",shift);
  
  /* The window */
  for(i = shift, j=0; j < window->size; j++,i++) {
    x = (i-window->size/2)/window->size;
    window->Y[j] =  0.42+0.5*cos(2*M_PI*x)+0.08*cos(4*M_PI*x);
    window->Y[j] *= window->Y[j];
    energy += window->Y[j]*window->Y[j];
  }
  
  /* Putting the first point to zero */
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  factor = 1/sqrt(energy);
  for(j=0;j<window->size;j++) {
    window->Y[j] *= factor;
  }
  return(factor);
}
/*************************/
/*    	Exponential      */
/*************************/

/*  Exponential asymmetric window */

static LWFLOAT Exponential(SIGNAL window,LWFLOAT shift)
{
  LWFLOAT decay = 1e4;
  int j;
  LWFLOAT i;
  LWFLOAT energy = 0.0;
  LWFLOAT factor;
  LWFLOAT a;
  LWFLOAT expon;
  
  /* Checking argument */
  if(window == NULL)
    Errorf("Asymmetric : NULL window");
  if(shift < 0 || shift >= 1)
    Errorf("Asymmetric : bad shift %g",shift);
  
  /* Damping factor */
  
  a = log(decay);
  expon = a/window->size;  /* scaled */
  
  
  /* The window */
  
  for(j=0, i=-1; i < window->size-1; j++, i++) {
    window->Y[j] = exp(-expon*i);
    energy += window->Y[j]*window->Y[j];
  }
  
  
  /* Putting the first point to zero*/
  
  energy -= window->Y[0]*window->Y[0];
  window->Y[0] = 0.0;
  
  /* Normalizing */
  
  factor = 1/sqrt(energy);
  
  for(j = 0; j < window->size; j++) {
    window->Y[j] *= factor;
  }
  /*  Printf("Asymmetric Window:(a: %f,factor: %f,size: %d)\n",a,factor,window->size); javi */
  
  return(factor);
}

/*************************/
/*          FoF	         */
/*************************/

/*  FoF (Fonction d'onde Formantique) window 			   */
/* 		g(t)= 0.5*exp(-a*t)*(1-cos(b*t))     0<= t <= PI/b */
/*		g(t)=     exp(-a*t)                      t >  PI/b */

LWFLOAT decayFoF = 1e5; /* 1e5 */ 
LWFLOAT betaFoF = M_PI/0.25; /* pi/0.25 */

static LWFLOAT FoF(SIGNAL window,LWFLOAT shift)
{