/[marvin]/marvin/src/libsip/algo/mfcc.cc
ViewVC logotype

Diff of /marvin/src/libsip/algo/mfcc.cc

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.4 by exos, Thu Sep 11 22:45:21 2003 UTC revision 1.5 by exos, Fri Sep 12 01:50:01 2003 UTC
# Line 39  Line 39 
39  *  *
40  *******************************************************************/  *******************************************************************/
41  #include "mfcc.hh"  #include "mfcc.hh"
42    #include "fft.hh"
43    
44  mfcc::mfcc ( ContSignal signal ){  mfcc::mfcc ( ContSignal signal ){
45  }  }
46    
47  mfcc::mfcc ( DiscSignal signal )  static int my_rint(double x)
48  {  {
49      if ( 2*x == (double)rint(2*x))
50        x+=0.0001;
51    
52      return((int)rint(x));
53    }
54    
55    mfcc::mfcc ( DiscSignal &signal )
56    {
57      //param: taille init, taille windowing
58    
59      double
60        lowestFrequency = 133.3333,
61        linearFilters = 13,
62        linearSpacing = 66.66666666,
63        logFilters = 27,
64        logSpacing = 1.0711703,
65        fftSize = 512,
66        cepstralCoefficients = 13,
67        windowSize = 256,
68        samplingRate = 1025,
69        totalFilters,
70        *mfccFilterWeights,
71        *freqs,
72        *triangleHeight,
73        *lower,
74        *upper,
75        *center,
76        *fftFreqs,
77        *hamWindow,
78        *mfccDCTMatrix,
79        *preEmphasized,
80        first,
81        last,
82        *fftData,
83        *empreinte,
84        maxi,
85        *earMag, *ceps, loga, windowStep, cols, frameRate =
86        100;
87      
88      int   i, j, k;
89      
90      preEmphasized = (double *)calloc(1024, sizeof(double));
91      totalFilters = linearFilters + logFilters;
92      freqs = (double *)calloc(sizeof(double), (int)totalFilters + 2);
93      for (i = 0; i < linearFilters; i++)
94        freqs[i] = lowestFrequency + i * linearSpacing;
95      for (i = (int)linearFilters; i < (int)totalFilters + 2; i++)
96        freqs[i] =
97          freqs[(int) linearFilters - 1] * pow(logSpacing,i - linearFilters + 1);
98      lower = (double*)calloc((int)totalFilters, sizeof(double));
99      upper = (double*)calloc((int)totalFilters, sizeof(double));
100      center = (double*)calloc((int)totalFilters, sizeof(double));
101      triangleHeight = (double*)calloc((int)totalFilters, sizeof(double));
102      fftData = (double*)calloc(sizeof(double), (int)fftSize);
103    
104      for (i = 0; i < (int)totalFilters; i++)
105        lower[i] = freqs[i];
106      for (i = 1; i < (int)totalFilters + 1; i++)
107        center[i - 1] = freqs[i];
108      for (i = 2; i < (int)totalFilters + 2; i++)
109        upper[i - 2] = freqs[i];
110    
111      hamWindow = (double*)calloc((int)windowSize, sizeof(double));
112      mfccFilterWeights = (double *)calloc((int)(totalFilters * fftSize), sizeof(double));
113      for (i = 0; i < totalFilters; i++)
114        triangleHeight[i] = 2 / (upper[i] - lower[i]);
115      fftFreqs = (double *)calloc((int)fftSize, sizeof(double));
116      for (i = 0; i < fftSize; i++)
117        fftFreqs[i] = (i / fftSize) * samplingRate;
118      for (i = 0; i < totalFilters; i++)
119        for (j = 0; j < (int)fftSize; j++) {
120          mfccFilterWeights[(int) fftSize * i + j] =
121            ((fftFreqs[j] > lower[i]) && (fftFreqs[j] <= center[i]))
122            * triangleHeight[i] * (fftFreqs[j] - lower[i]) / (center[i] - lower[i]) +
123            ((fftFreqs[j] > center[i]) && (fftFreqs[j] <upper[i])) * triangleHeight[i] * (upper[i] - fftFreqs[j]) /
124            (upper[i] - center[i]);
125        }
126      //    for (i=0;i<40*512;i++) printf("%d %d %f\n",i/512,i%512,mfccFilterWeights[i]);
127      for (i = 0; i < windowSize; i++)
128        hamWindow[i] = 0.54 - 0.46 * cos(2 * PI * i / windowSize);
129      mfccDCTMatrix = (double*)
130        calloc(sizeof(double), (int)(cepstralCoefficients * totalFilters));
131      for (i = 0; i < cepstralCoefficients; i++)
132        for (j = 0; j < totalFilters; j++)
133          mfccDCTMatrix[j * (int) cepstralCoefficients + i] =
134            1 / sqrt(totalFilters / 2) * cos(i * (2 * j + 1) * PI / 2 / totalFilters);
135      //    for (j=0;j<13*40;j++) printf("%d %d %f \n",j/13,j%13,mfccDCTMatrix[j]);
136      for (i = 0; i < totalFilters; i++)
137        mfccDCTMatrix[i * (int) cepstralCoefficients] *= 1 / sqrt(2);
138      for (i = 1; i < 1024; i++)
139        preEmphasized[i] = signal[i] + signal[i - 1] * -0.97;
140      preEmphasized[0] = signal[0];
141      //for (i=0;i<1024;i++) printf("%d %f \n",i,preEmphasized[i]);
142      earMag = (double*)malloc(sizeof(double) * (int)totalFilters);
143      windowStep = samplingRate / frameRate;
144      cols = (int) ((1024 - windowSize) / windowStep);      
145      //printf("%d",(int)cols);
146      ceps = (double*)malloc((int)(sizeof(double) * cols * linearFilters));
147      for (i = 0; i < cols; i++) {
148        first = i * windowStep + 1;
149        last = first + windowSize - 1;
150        for (j = 0; j < fftSize; j++)
151          fftData[j] = 0;
152        for (j = 0; j < windowSize; j++)
153          fftData[j] =
154            preEmphasized[(int) my_rint(first + j - 1)] *
155            hamWindow[j];
156        //for (j=0; j<512;j++) printf("%d %f\n",j,fftData[j]);
157        /*    for (j = 0; j < fftSize; j++) {
158              inR[j] = fftData[j];
159              inI[j] = 0;
160              }
161              fft_double(512, 0, inR, inI, ouR, ouI);
162              for (j = 0; j < fftSize; j++)
163              fftMag[j] = sqrt(ouR[j] * ouR[j] + ouI[j] * ouI[j]);
164              //for (j=0; j<512;j++) printf("%d %f\n",j,fftMag[j]);
165              */
166        DiscSignal fftMag = (new fft(*new DiscSignal(fftData, (int)fftSize), true))->get_magnitude();
167        for (j = 0; j < totalFilters; j++) {
168          loga = 0;
169          for (k = 0; k < fftSize; k++)
170            loga +=
171              fftMag[k] * mfccFilterWeights[k + j * (int) fftSize];
172          earMag[j] = log10(loga);
173        }
174        //              for (j=0;j<40;j++) printf("%d %f\n",j,earMag[j]);
175        for (j = 0; j < linearFilters; j++) {
176          loga = 0;
177          for (k = 0; k < totalFilters; k++)
178            loga +=
179              earMag[k] * mfccDCTMatrix[j + k * (int) linearFilters];
180          ceps[j + i * (int) linearFilters] = loga;
181        }
182      }
183      //////////////////////moyenne//////////////////////
184      empreinte = (double*)malloc(sizeof(double) * 12);
185      //    for (i=0;i<68;i++)
186      //      printf("%d %d %f\n",i/13,i%13,ceps[i]);
187      //            for (i =0; i<13;i++) for(j=0;j<6;j++) printf("%d %d %f\n",j,i,ceps[i + j*13]);
188      for (i = 1; i < 13; i++)
189        empreinte[i - 1] = (ceps[i] + ceps[i + 13] + ceps[i+26] + ceps[i+39] + ceps[i+52] + ceps[i+65]);
190    
191      maxi=0;
192      for (i=0;i<12;i++) if (maxi*maxi<empreinte[i]*empreinte[i]) maxi=empreinte[i];
193      for (i=0;i<12;i++) empreinte[i]=(empreinte[i]/fabs(maxi))*2;
194      (new DiscSignal(empreinte, 13))->perl_capture();
195      //for (i=0;i<12;i++) printf("%d %lf\n",i,empreinte[i]);
196      free(mfccFilterWeights);
197      free(triangleHeight);
198      free(lower);
199      free(upper);
200      free(center);
201      free(fftFreqs);
202      free(hamWindow);
203      free(mfccDCTMatrix);
204      free(preEmphasized);
205      free(fftData);
206    
207      free(earMag);
208      free(ceps);
209    
210    
211    
212      //    return (empreinte);
213    
214    
215  }  }

Legend:
Removed from v.1.4  
changed lines
  Added in v.1.5

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26