libghost/quant_sinusoids.c

   1 /**
   2    @file sinusoids.c
   3    @brief Sinusoid extraction/synthesis
   4 */
   5
   6 /* Copyright (C) 2005
   7
   8    This library is free software; you can redistribute it and/or
   9    modify it under the terms of the GNU Lesser General Public
  10    License as published by the Free Software Foundation; either
  11    version 2.1 of the License, or (at your option) any later version.
  12
  13    This library is distributed in the hope that it will be useful,
  14    but WITHOUT ANY WARRANTY; without even the implied warranty of
  15    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  16    Lesser General Public License for more details.
  17
  18    You should have received a copy of the GNU Lesser General Public
  19    License along with this library; if not, write to the Free Software
  20    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  21 */
  22
  23
  24 #include <math.h>
  25 #include "sinusoids.h"
  26 #include <stdio.h>
  27 #include "ghost.h"
  28
  29 #define MIN(a,b) ((a)<(b) ? (a):(b))
  30 #define MAX(a,b) ((a)>(b) ? (a):(b))
  31
  32
  33
  34
  35 /* Quantizes the extracted sinusoids
  36  * w are the frequencies in current frame
  37  * N is the number of sinusoids
  38  * SineStruct is structure holding quantized sine params
  39 */
  40
  41 void quantize_sinusoids(float *wi,int N,QuantSine SineStruct[SINUSOIDS])
  42 {
  43
  44 /* sort the current sinusoids */
  45         float w[SINUSOIDS];
  46         int i,j,lindex = 0;
  47         static QuantSine prevSineStruct[SINUSOIDS];
  48
  49
  50         for(i=0;i<N;i++)
  51         {
  52                 w[i] = wi[i];
  53         }
  54
  55         for(i=0;i<N;i++)
  56         {
  57                 for(j=(i+1);j<N;j++)
  58                 {
  59                         if(w[j]< w[i])
  60                         {
  61                                 float t = w[i];
  62                                 w[i] = w[j];
  63                                 w[j] = t;
  64                         }
  65                 }
  66         }
  67
  68         /* start with lowest freq and group its harmonics. Go on till
  69 all sinuosids belong to a group*/
  70
  71         for(i=0;i<N;i++)
  72         {
  73                 if(w[i]!= -1)
  74                 {
  75                         SineStruct[lindex].basefreq = w[i];
  76                         SineStruct[lindex].harmonicstate;
  77                         for(j=(i+1);j<N;j++)
  78                         {
  79                                 if(w[j]!= -1)
  80                                 {
  81                                         int harmonic = w[j]/SineStruct[lindex].basefreq;
  82
  83                                         if((w[j]/SineStruct[lindex].basefreq - (float)(harmonic) < 0.01)&&
  84                                                 (w[j]/SineStruct[lindex].basefreq - (float)(harmonic) > -0.01))
  85                                         {
  86                                                 w[j] = -1;
  87                                                 SineStruct[lindex].harmonicstate = (SineStruct[lindex].harmonicstate)
  88                                                                                                                         | (1<<harmonic);
  89                                         }
  90
  91                                 }
  92                         }
  93                         lindex++;
  94                 }
  95         }
  96
  97
  98 /* compare the current groups with groups of prev frames
  99 and see which all have to be transmitted */
 100         i= 0;
 101         while(i<lindex)
 102         {
 103                 /* compare the current set with all the freq in the prev bands*/
 104                 if(SineStruct[i].basefreq < prevSineStruct[j].basefreq)
 105                 {
 106                         //add i
 107                         i++;
 108                 }
 109                 else if(SineStruct[i].basefreq > prevSineStruct[j].basefreq)
 110                 {
 111                         //rem j
 112                         j++;
 113                 }
 114                 else
 115                 {
 116                         if(SineStruct[i].harmonicstate != prevSineStruct[j].harmonicstate)//state[i] != prevstate [j]
 117                         {
 118                                 //update state of i
 119                         }
 120                 }
 121
 122         }
 123 }
 124
 125 /* Inverse quantizes the extracted sinusoids
 126  * w are the frequencies in current frame
 127  * N is the number of sinusoids
 128  * len is the frame size
 129 */
 130
 131 void invquantize_sinusoids(float *x,float *w,int N,int len)
 132 {
 133
 134 /*add the sinuosids newly sent in current frame to wi vector */
 135
 136         /* check for all teh sinusoids; add the ones which are to be added and remove the ones which are to be removed to the
 137         w vector already existing */
 138
 139 }
 140
 141 /* Synthesizes the signal from the sinusoids
 142  * w are the frequencies
 143  * ai are the cos(x) coefficients
 144  * bi are the sin(x) coefficients
 145  * y is the synthesized signal by summing all the params
 146  * N is the number of sinusoids
 147  * len is the frame size
 148 */
 149 void generate_sinusoids(float *w,
 150                                                 float *window,
 151                                                 float *ai,
 152                                                 float *bi,
 153                                                 float *y,
 154                                                 int N,
 155                                                 int len)
 156 {
 157         int i,j;
 158         float cos_table[SINUSOIDS][LENGTH];
 159     float sin_table[SINUSOIDS][LENGTH];
 160
 161         for (i=0;i<N;i++)
 162     {
 163       float tmp1=0, tmp2=0;
 164       float tmp3=0, tmp4=0;
 165       for (j=0;j<len;j++)
 166       {
 167          cos_table[i][j] = cos(w[i]*j)*window[j];
 168          sin_table[i][j] = sin(w[i]*j)*window[j];
 169       }
 170     }
 171
 172         for (i=0;i<N;i++)
 173     {
 174          for (j=0;j<len;j++)
 175          {
 176             y[j] += ai[i]*cos_table[i][j] + bi[i]*sin_table[i][j];
 177          }
 178     }
 179
 180 }