mp_graph.c

LastWave

    AppendFloat2Listv(lv,o->rh);
    return(YES);
  }

  return(NO);
}


/*********************************************************************
 * Functions to fill the display area with the Wigner distribution
 *********************************************************************/


/*
 * The Wigner-Ville distribution of a Gaussian atom is
 *  GWV(2*(u/2^o)) x GWV(2*(2*pi*\sigma^2* k*2^o/GABOR_MAX_FREQID))
 *
 *  GWV(x) = e^{-x^2/4\sigma^2}
 */
LWFLOAT *GaussianWignerVille(int sizeGWV)
{
    LWFLOAT *GWV;
    LWFLOAT c = 1/(4*theGaussianSigma2*sizeGWV*sizeGWV);
    int i;
        
    /* Allocation */
    if ((GWV = FloatAlloc(sizeGWV)) == NULL) 
        Errorf("GaussianWignerVilleTime : Mem. Alloc. failed!");
        
    /* Computation of e^{-x^2/\sigma^2}, x = i/sizeGWV*/
    for (i = 0; i < sizeGWV; i++) 
    {
        GWV[i] = exp(-c*i*i);
    }
    return(GWV);
}

/* 
 * The Wigner-Ville distribution of a FoF atom is
 *  FWV(2*(u/2^o)) x GWV(2*(2*pi*\sigma^2* k*2^o/GABOR_MAX_FREQID))
 *
 *  FWV(x) = ???
 */
LWFLOAT *FoFWignerVille(int sizeFWV)
{
    LWFLOAT *FWV;
    LWFLOAT a,expon,beta,max;
    int i;
        
    /* Allocation */
    if ((FWV = FloatAlloc(sizeFWV)) == NULL) 
        Errorf("FoFWignerVilleTime : Mem. Alloc. failed!");
    
    /* Damping factor */
    a = log(decayFoF);
    
    /* scale */
    expon = a/sizeFWV; 
    beta = betaFoF/sizeFWV;        
    
    max=0.0;

    /* Computation of FoF window */
    for (i = 0; i <= M_PI/beta; i++) 
    {
        FWV[i] = 0.5*(1-cos(beta*i))*exp(-expon*i);
	if (FWV[i] > max) max=FWV[i];
    }
    for (; i<sizeFWV; i++)
    {
	FWV[i] = exp(-expon*i);
	if (FWV[i] > max) max=FWV[i];
    }

    /* Normalizing to the maximum */
    for (i = 0; i <sizeFWV; i++) 
    {
        FWV[i]=FWV[i]/max;
    }
    return(FWV);
}

/* 
 * Add the Wigner-Ville time-frequency distribution
 * of an atom to a time-frequency 'map' of nRow*nCol
 */
void DisplayGaborAtom(GOBJECT obj,LWFLOAT *map,ATOM atom,int atomK,int mapx,int mapy,int nRow,int nCol)
{
  GRAPHBOOK graph;
  BOOK book;

  /* The templates of the Wigner-Ville distributions */
  static int sizeG = 256;
  static LWFLOAT *GWV = NULL;
  static LWFLOAT *FWV = NULL;

  LWFLOAT timeMin,timeWidth,freqMin,freqWidth;
  LWFLOAT nTimeIdPerCol,nFreqIdPerRow;
  
  char flagAsym;

  int x,y,xcenter,ycenter,dx,dy;
  int xcurrent;
  int xmin,xmax,ymin,ymax;
  
  LWFLOAT iWV;
  int  t,f,k;

  LWFLOAT timeWV,freqWV;
  

  /* Some inits */
  graph = (GRAPHBOOK) obj;
  book = graph->book;
  if (book == NULL) return;
  CheckBookNotEmpty(book);
  flagAsym = GetFlagAsymWindowShape(atom->windowShape);

  /* Is this atom to be displayed in this graph ? */
  if(!IsVisible(graph,atom)) return;

  /* Initialize the templates of Wigner-Ville distribution */
  if(flagAsym) {
    if(FWV == NULL) {
      FWV = FoFWignerVille(sizeG);
    }
    if(GWV == NULL) {
      GWV = GaussianWignerVille(sizeG);
    }
  }
  else {
    if(GWV == NULL) {
      GWV = GaussianWignerVille(sizeG);
    }
  }
 
  /* 
   * Get the pixel coordinates (xcenter,ycenter)
   * of the 'center' of the atom, in map[] :
   */
  /* the coordinates in the gobject */
  Local2Global(obj,TimeId2Time(book,atom->timeId),FreqId2Freq(book,atom->freqId),&xcenter,&ycenter);
  /* change the origin to the corner of the map */
  xcenter -=mapx;
  ycenter -=mapy;

  /*
   * Get the pixel rectangle (xmin,xmax,ymin,ymax) corresponding to the atom
   * as well as the conversion ratios between pixel and (time/freq)Id scales
   */
  /* (time/freq)(Min/Width) <=> size of the map in (time/freq) units */
  Global2LocalRect(obj,mapx,mapy,nCol,nRow,
		   &timeMin,&freqMin,&timeWidth,&freqWidth,NormalRect);
  /* Conversion ratios */
  nTimeIdPerCol = (Time2TimeId(book,timeMin+timeWidth)-Time2TimeId(book,timeMin))/nCol;
  nFreqIdPerRow = (Freq2FreqId(book,freqMin+freqWidth)-Freq2FreqId(book,freqMin))/nRow;

  /* The rectangle for a Gabor atom */
  /* dx and dy are HALF the width and heighth of the box (in pixels) */
  if (!flagAsym) {
    dx = (int) (atom->windowSize/(2*nTimeIdPerCol));
  }
  // Except for asymmetric atoms !
  else {
    dx = (int) (atom->windowSize/nTimeIdPerCol);
  }
  dy = (int) ((GABOR_MAX_FREQID/(4*M_PI*theGaussianSigma2*atom->windowSize))/nFreqIdPerRow);
  /* Now the box itself */
  xmax = xcenter+dx;
  if (!flagAsym) {
    xmin = xcenter-dx;
  }
  else {
    xmin = xcenter;
  }
  ymax = ycenter+dy;
  ymin = ycenter-dy;
  
  /* Is there anything to be displayed ?
   * i.e. does the box intersect the pixmap ? */
  if (xmax < 0 || ymax < 0 || xmin >= nCol || ymin >= nRow) return;

  /* We do not print outside the pixmap */
  xmin = MAX(0,xmin);
  xmax = MIN(nCol-1,xmax);
  ymin = MAX(0,ymin);
  ymax = MIN(nRow-1,ymax);

  /* 
   * Fill the pixmap
   */
  /* Case of a Gabor atom */
  /* columns x = 'time' */
  for (x = xmin; x <= xmax; x++) {
    /* 
     * Compute the index in the template of Wigner-Ville distrib. 
     */
    /* current distance in pixels from the 'time-center' of the atom */
    iWV = abs(x-xcenter);
    /* conversion in timeId units */
    iWV *= nTimeIdPerCol; 

    /* conversion in percentage of the time-support of the template WV */
    /* For symetric atoms, we have two equal sides of size 2^o/2 */
    if (!flagAsym) {
      iWV /= atom->windowSize/2;
    }
    /* For asymetric ones, only one side of size 2^o */
    else {
      iWV /= atom->windowSize;
    }

    /* conversion in index in the template WV */
    t = (int) (sizeG*iWV);
    
    /* indexes outside the tabulated template WV 
     * are expected to be be very small and are not displayed
     */
    if (t>=sizeG) continue;
    
    /* 
     * It is time to use the template WV to get
     * the amplitude of the time factor of the WV.
     */
    if(flagAsym) {
      timeWV = FWV[t];
    }
    else {
      timeWV = GWV[t];
    }

    /* Now, we take into account the CHIRP to compute
     * (in pixel coordinates) the instantaneous frequency 
     * f = atomtime+atomchirp*(delta t)
     * with the 'delta t' corresponding to the current value of 'x'
     */
    Local2Global(obj,
                 TimeId2Time(book,atom->timeId+(x-xcenter)*nTimeIdPerCol),
                 FreqId2Freq(book,atom->freqId+atom->chirpId*(x-xcenter)*nTimeIdPerCol),&xcurrent,&ycenter); 
    xcurrent -=mapx;
    ycenter  -=mapy;  
    
    /* compute the local frequency box */
    ymax = ycenter+dy;
    ymin = ycenter-dy;
    ymin = MAX(0,ymin);
    ymax = MIN(nRow-1,ymax);
        
    /* loop on frequency/row within the local box */
    for (y = ymin; y <= ymax; y++) {
      /* 
       * Compute the index in the template of Wigner-Ville distrib. 
       */
      /* current distance in pixels from the 'frequency-center' of the atom */
      iWV = abs(y-ycenter);
      /* conversion in freqId units */
      iWV *= nFreqIdPerRow; 
      /* conversion in percentage of the freq-support of the template WV */
      iWV *= (4*M_PI*theGaussianSigma2*atom->windowSize)/(GABOR_MAX_FREQID);
      /* conversion in index in the template WV */
      f = (int) (sizeG*iWV);
            
      /* indexes outside the tabulated template WV 
       * are expected to be be very small and are not displayed
       */
      if (f>=sizeG) continue;
            
      /* 
       * It is time to use the template WV to get
       * the amplitude of the freq factor of the WV.
       */
      freqWV = GWV[f];
      /* Adding contribution at the right position in 'map' */
      k = y*nCol+x;
      map[k] += atom->coeff2*timeWV*freqWV;
    }
  }
  return;
}



/* The drawing procedure */
static void _DrawGraphBook (WINDOW win, GOBJECT obj, int x, int y,int w,int h)
{
  GRAPHBOOK graph;
  BOOK book;

  unsigned short kMin,kMax;
  int nColors;

  int nRow,nCol;
  int color;

  static LWFLOAT *map         = NULL;
  static int   mapAllocSize = 0;

  LWFLOAT maxMap,maxPartialEnergy,moleculeEnergy;

  unsigned long n;
  MOLECULE molecule;
  
  ATOM atom;

  int i,j;
  unsigned long k;


  /* Some inits */
  graph = (GRAPHBOOK) obj;
  book  = graph->book;
  if (book == NULL) return;
  CheckBookNotEmpty(book);

  /* Get the number of colors of the colormap */
  nColors = ColorMapSize(graph->cm);

  /* Allocation */
  nCol = w;
  nRow = h;
  WInitPixMap(nRow,nCol);
  
  /* Allocate a 'local map' of amplitudes and initialize with zeros */
  if(map != NULL && mapAllocSize < nRow*nCol) {
    Free(map);
    map          = NULL;
    mapAllocSize = 0;
  }
  if(map == NULL) {
    if ((map = FloatAlloc(nRow*nCol)) == NULL) 
        Errorf("_DrawGraphBook : Mem. Alloc. failed!");
    mapAllocSize = nRow*nCol;
  }
  for(k = 0; k < mapAllocSize; k++) map[k] = 0;

  /*
   * Do the plot in the local map
   */
  maxMap = -1;

  /* Draw each molecule within the graph's limits */
  for(n = graph->nMin; n <= MIN(book->size-1,graph->nMax); n++)  {
    molecule = GetBookMolecule(book,n);
    switch(graph->flagFund) {
      /* Case when we display all partials within a certain range */
    case FundAll :
      kMin = 0;
      kMax = molecule->dim-1;
      break;
      /* Case when we display only the most energetic partial */
    case FundMax :
    case FundSum :
      maxPartialEnergy = -1;
      moleculeEnergy       = 0;    
      /*
       * We locate the most energetic partial and compute
       * maxPartialEnergy    = the energy of this partial
       * moleculeEnergy          = the total energy of the molecule
       */
      for(k = 0; k < molecule->dim; k++) {
        atom = GetMoleculeAtom(molecule,0,k);
        moleculeEnergy += atom->coeff2;
        if (maxPartialEnergy < atom->coeff2) {
          maxPartialEnergy = MAX(maxPartialEnergy,atom->coeff2);
          kMin = kMax = k;
        }
      }
      break;
    }
      
    /* Now we do display the desired partials */
    for(k = kMin; k <= kMax; k++) {
      atom = GetMoleculeAtom(molecule,0,k);
      /*
       * Case 'FundSum' : 
       * the most energetic with the total energy of the molecule 
       */
      if (graph->flagFund == FundSum) {
	/* We temporarily modify the coeff2 and memorize the right one
           in maxPartialEnergy */
        maxPartialEnergy = atom->coeff2;
        atom->coeff2 = moleculeEnergy;
      }
             
      maxMap = MAX (maxMap,atom->coeff2);
          
      DisplayGaborAtom(obj,map,atom,k,x,y,nRow,nCol);
      
      /* We restore atom->coeff2 if necessary */
      if (graph->flagFund == FundSum) atom->coeff2 = maxPartialEnergy;

    }
  }
  
  /* Normalize and replace the map by its log (dB) or power, pixel by pixel */
  if(graph->flagDb) {
    for(k = 0; k < nCol*nRow; k++) {
      if(map[k] != 0.0) {
	/* Normalize */
	map[k] /= maxMap;
	/* Take decibels */
	map[k] = 10*log10(map[k]);
	/* Scale the dynamics with 'exponent' */
	map[k] /= graph->exponent;
	map[k] += 1;
	map[k] = MAX(map[k],0);
      }
    }
  }    
  else {
    for(k = 0; k < nCol*nRow; k++) { 
       if(map[k] != 0.0) {
	/* Normalize */
	map[k] /= maxMap;
	map[k] = pow(map[k],graph->exponent);
      }
    }
  }
  
  /* Puts the 'local map' into the pixmap, using the color-code */
  /* rows = 'freq' */
  for(i = 0; i < nRow; i++) {
    /* columns =  'time' */
    for(j = 0; j < nCol; j++) {
      k = i*nCol+j;
      color = (int) (nColors * map[k]);
      color = (color>=nColors ? nColors-1 : color);
      color = (color<0 ? 0 : color);
      color += graph->cm;
          
      WSetPixelPixMap(i,j,color);
    }
  }
  WDisplayPixMap(win,x,y);   
}

/* Defining the GraphBook gclass */  
void DefineGraphBook(void)
{
  theGraphBookClass          = NewGClass("GraphBook",theGObjectClass,"mp"); 
  theGraphBookClass->nbBytes = sizeof(GraphBook);
  theGraphBookClass->init    = _InitGraphBook;
  theGraphBookClass->deleteContent = _DeleteContentGraphBook;
  theGraphBookClass->set     = _SetGraphBook;
  theGraphBookClass->draw    = _DrawGraphBook;   
  theGraphBookClass->varType = bookType;
  theGraphBookClass->flags  &= ~(GClassMoveResize+GClassLocalCoor);
  theGraphBookClass->info    = "Graphic Class that allows to display Book";
}


/* EOF */