diff options
| -rw-r--r-- | sms/cepstrum.c | 278 |
1 files changed, 139 insertions, 139 deletions
diff --git a/sms/cepstrum.c b/sms/cepstrum.c index 26af068..8a56787 100644 --- a/sms/cepstrum.c +++ b/sms/cepstrum.c @@ -33,45 +33,45 @@ typedef struct { - int nPoints; - int nCoeff; - gsl_matrix *pM; - gsl_matrix *pMt; - gsl_matrix *pR; - gsl_matrix *pMtMR; - gsl_vector *pXk; - gsl_vector *pMtXk; - gsl_vector *pC; - gsl_permutation *pPerm; + int nPoints; + int nCoeff; + gsl_matrix *pM; + gsl_matrix *pMt; + gsl_matrix *pR; + gsl_matrix *pMtMR; + gsl_vector *pXk; + gsl_vector *pMtXk; + gsl_vector *pC; + gsl_permutation *pPerm; } CepstrumMatrices; void FreeDCepstrum(CepstrumMatrices *m) { - gsl_matrix_free(m->pM); - gsl_matrix_free(m->pMt); - gsl_matrix_free(m->pR); - gsl_matrix_free(m->pMtMR); - gsl_vector_free(m->pXk); - gsl_vector_free(m->pMtXk); - gsl_vector_free(m->pC); - gsl_permutation_free (m->pPerm); + gsl_matrix_free(m->pM); + gsl_matrix_free(m->pMt); + gsl_matrix_free(m->pR); + gsl_matrix_free(m->pMtMR); + gsl_vector_free(m->pXk); + gsl_vector_free(m->pMtXk); + gsl_vector_free(m->pC); + gsl_permutation_free (m->pPerm); } void AllocateDCepstrum(int nPoints, int nCoeff, CepstrumMatrices *m) { - if(m->nPoints != 0 || m->nCoeff != 0) - FreeDCepstrum(m); - m->nPoints = nPoints; - m->nCoeff = nCoeff; - m->pM = gsl_matrix_alloc(nPoints, nCoeff); - m->pMt = gsl_matrix_alloc(nCoeff, nPoints); - m->pR = gsl_matrix_calloc(nCoeff, nCoeff); - m->pMtMR = gsl_matrix_alloc(nCoeff, nCoeff); - m->pXk = gsl_vector_alloc(nPoints); - m->pMtXk = gsl_vector_alloc(nCoeff); - m->pC = gsl_vector_alloc(nCoeff); - m->pPerm = gsl_permutation_alloc (nCoeff); + if(m->nPoints != 0 || m->nCoeff != 0) + FreeDCepstrum(m); + m->nPoints = nPoints; + m->nCoeff = nCoeff; + m->pM = gsl_matrix_alloc(nPoints, nCoeff); + m->pMt = gsl_matrix_alloc(nCoeff, nPoints); + m->pR = gsl_matrix_calloc(nCoeff, nCoeff); + m->pMtMR = gsl_matrix_alloc(nCoeff, nCoeff); + m->pXk = gsl_vector_alloc(nPoints); + m->pMtXk = gsl_vector_alloc(nCoeff); + m->pC = gsl_vector_alloc(nCoeff); + m->pPerm = gsl_permutation_alloc (nCoeff); } /*! \brief Discrete Cepstrum Transform @@ -98,58 +98,58 @@ void AllocateDCepstrum(int nPoints, int nCoeff, CepstrumMatrices *m) * \param iMaxFreq maximum frequency of cepstrum */ void sms_dCepstrum( int sizeCepstrum, sfloat *pCepstrum, int sizeFreq, sfloat *pFreq, sfloat *pMag, - sfloat fLambda, int iMaxFreq) + sfloat fLambda, int iMaxFreq) { - int i, k; - sfloat factor; - sfloat fNorm = PI / (sfloat)iMaxFreq; /* value to normalize frequencies to 0:0.5 */ - //static sizeCepstrumStatic - static CepstrumMatrices m; - //printf("nPoints: %d, nCoeff: %d \n", m.nPoints, m.nCoeff); - if(m.nPoints != sizeCepstrum || m.nCoeff != sizeFreq) - AllocateDCepstrum(sizeFreq, sizeCepstrum, &m); - int s; /* signum: "(-1)^n, where n is the number of interchanges in the permutation." */ - /* compute matrix M (eq. 4)*/ - for (i=0; i<sizeFreq; i++) - { - gsl_matrix_set (m.pM, i, 0, 1.); // first colum is all 1 - for (k=1; k <sizeCepstrum; k++) - gsl_matrix_set (m.pM, i, k , 2.*sms_sine(PI_2 + fNorm * k * pFreq[i]) ); - } + int i, k; + sfloat factor; + sfloat fNorm = PI / (sfloat)iMaxFreq; /* value to normalize frequencies to 0:0.5 */ + //static sizeCepstrumStatic + static CepstrumMatrices m; + //printf("nPoints: %d, nCoeff: %d \n", m.nPoints, m.nCoeff); + if(m.nPoints != sizeCepstrum || m.nCoeff != sizeFreq) + AllocateDCepstrum(sizeFreq, sizeCepstrum, &m); + int s; /* signum: "(-1)^n, where n is the number of interchanges in the permutation." */ + /* compute matrix M (eq. 4)*/ + for (i=0; i<sizeFreq; i++) + { + gsl_matrix_set (m.pM, i, 0, 1.); // first colum is all 1 + for (k=1; k <sizeCepstrum; k++) + gsl_matrix_set (m.pM, i, k , 2.*sms_sine(PI_2 + fNorm * k * pFreq[i]) ); + } - /* compute transpose of M */ - gsl_matrix_transpose_memcpy (m.pMt, m.pM); - - /* compute R diagonal matrix (for eq. 7)*/ - factor = COEF * (fLambda / (1.-fLambda)); /* \todo why is this divided like this again? */ - for (k=0; k<sizeCepstrum; k++) - gsl_matrix_set(m.pR, k, k, factor * powf((sfloat) k,2.)); + /* compute transpose of M */ + gsl_matrix_transpose_memcpy (m.pMt, m.pM); - /* MtM = Mt * M, later will add R */ - gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., m.pMt, m.pM, 0.0, m.pMtMR); - /* add R to make MtMR */ - gsl_matrix_add (m.pMtMR, m.pR); + /* compute R diagonal matrix (for eq. 7)*/ + factor = COEF * (fLambda / (1.-fLambda)); /* \todo why is this divided like this again? */ + for (k=0; k<sizeCepstrum; k++) + gsl_matrix_set(m.pR, k, k, factor * powf((sfloat) k,2.)); - /* set pMag in X and multiply with Mt to get pMtXk */ - for(k = 0; k <sizeFreq; k++) - gsl_vector_set(m.pXk, k, log(pMag[k])); - gsl_blas_dgemv (CblasNoTrans, 1., m.pMt, m.pXk, 0., m.pMtXk); + /* MtM = Mt * M, later will add R */ + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., m.pMt, m.pM, 0.0, m.pMtMR); + /* add R to make MtMR */ + gsl_matrix_add (m.pMtMR, m.pR); - /* solve x (the cepstrum) in Ax = b, where A=MtMR and b=pMtXk */ + /* set pMag in X and multiply with Mt to get pMtXk */ + for(k = 0; k <sizeFreq; k++) + gsl_vector_set(m.pXk, k, log(pMag[k])); + gsl_blas_dgemv (CblasNoTrans, 1., m.pMt, m.pXk, 0., m.pMtXk); - /* ==== the Cholesky Decomposition way ==== */ - /* MtM is 'symmetric and positive definite?' */ - //gsl_linalg_cholesky_decomp (m.pMtMR); - //gsl_linalg_cholesky_solve (m.pMtMR, m.pMtXk, m.pC); + /* solve x (the cepstrum) in Ax = b, where A=MtMR and b=pMtXk */ - /* ==== the LU decomposition way ==== */ - gsl_linalg_LU_decomp (m.pMtMR, m.pPerm, &s); - gsl_linalg_LU_solve (m.pMtMR, m.pPerm, m.pMtXk, m.pC); + /* ==== the Cholesky Decomposition way ==== */ + /* MtM is 'symmetric and positive definite?' */ + //gsl_linalg_cholesky_decomp (m.pMtMR); + //gsl_linalg_cholesky_solve (m.pMtMR, m.pMtXk, m.pC); - - /* copy pC to pCepstrum */ - for(i = 0; i < sizeCepstrum; i++) - pCepstrum[i] = gsl_vector_get (m.pC, i); + /* ==== the LU decomposition way ==== */ + gsl_linalg_LU_decomp (m.pMtMR, m.pPerm, &s); + gsl_linalg_LU_solve (m.pMtMR, m.pPerm, m.pMtXk, m.pC); + + + /* copy pC to pCepstrum */ + for(i = 0; i < sizeCepstrum; i++) + pCepstrum[i] = gsl_vector_get (m.pC, i); } /*! \brief Spectrum Envelope from Cepstrum @@ -163,36 +163,36 @@ void sms_dCepstrum( int sizeCepstrum, sfloat *pCepstrum, int sizeFreq, sfloat *p */ void sms_dCepstrumEnvelope(int sizeCepstrum, sfloat *pCepstrum, int sizeEnv, sfloat *pEnv) { - - static sfloat *pFftBuffer; - static int sizeFftArray = 0; - int sizeFft = sizeEnv << 1; - int i; + + static sfloat *pFftBuffer; + static int sizeFftArray = 0; + int sizeFft = sizeEnv << 1; + int i; + if(sizeFftArray != sizeFft) + { + if(sizeFftArray != 0) free(pFftBuffer); + sizeFftArray = sms_power2(sizeFft); if(sizeFftArray != sizeFft) { - if(sizeFftArray != 0) free(pFftBuffer); - sizeFftArray = sms_power2(sizeFft); - if(sizeFftArray != sizeFft) - { - sms_error("bad fft size, incremented to power of 2"); - } - if ((pFftBuffer = (sfloat *) malloc(sizeFftArray * sizeof(sfloat))) == NULL) - { - sms_error("could not allocate memory for fft array"); - return; - } + sms_error("bad fft size, incremented to power of 2"); + } + if ((pFftBuffer = (sfloat *) malloc(sizeFftArray * sizeof(sfloat))) == NULL) + { + sms_error("could not allocate memory for fft array"); + return; } - memset(pFftBuffer, 0, sizeFftArray * sizeof(sfloat)); + } + memset(pFftBuffer, 0, sizeFftArray * sizeof(sfloat)); - pFftBuffer[0] = pCepstrum[0] * 0.5; - for (i = 1; i < sizeCepstrum-1; i++) - pFftBuffer[i] = pCepstrum[i]; + pFftBuffer[0] = pCepstrum[0] * 0.5; + for (i = 1; i < sizeCepstrum-1; i++) + pFftBuffer[i] = pCepstrum[i]; - sms_fft(sizeFftArray, pFftBuffer); + sms_fft(sizeFftArray, pFftBuffer); - for (i = 0; i < sizeEnv; i++) - pEnv[i] = powf(EXP, 2. * pFftBuffer[i*2]); + for (i = 0; i < sizeEnv; i++) + pEnv[i] = powf(EXP, 2. * pFftBuffer[i*2]); } /*! \brief main function for computing spectral envelope from sinusoidal peaks @@ -206,54 +206,54 @@ void sms_dCepstrumEnvelope(int sizeCepstrum, sfloat *pCepstrum, int sizeEnv, sfl */ void sms_spectralEnvelope( SMS_Data *pSmsData, SMS_SEnvParams *pSpecEnvParams) { - int i, k; - int sizeCepstrum = pSpecEnvParams->iOrder+1; - //int nPeaks = 0; - static sfloat pFreqBuff[1000], pMagBuff[1000]; - - /* \todo see if this memset is even necessary, once working */ - //memset(pSmsData->pSpecEnv, 0, pSpecEnvParams->nCoeff * sizeof(sfloat)); + int i, k; + int sizeCepstrum = pSpecEnvParams->iOrder+1; + //int nPeaks = 0; + static sfloat pFreqBuff[1000], pMagBuff[1000]; - /* try to store cepstrum coefficients in pSmsData->nEnvCoeff always. - if cepstrum is what is wanted, memset the rest. otherwise, hand this array 2x to dCepstrumEnvelope */ - if(pSpecEnvParams->iOrder + 1> pSmsData->nEnvCoeff) - { - sms_error("cepstrum order is larger than the size of the spectral envelope"); - return; - } + /* \todo see if this memset is even necessary, once working */ + //memset(pSmsData->pSpecEnv, 0, pSpecEnvParams->nCoeff * sizeof(sfloat)); - /* find out how many tracks were actually found... many are zero - \todo is this necessary? */ - for(i = 0, k=0; i < pSmsData->nTracks; i++) + /* try to store cepstrum coefficients in pSmsData->nEnvCoeff always. + if cepstrum is what is wanted, memset the rest. otherwise, hand this array 2x to dCepstrumEnvelope */ + if(pSpecEnvParams->iOrder + 1> pSmsData->nEnvCoeff) + { + sms_error("cepstrum order is larger than the size of the spectral envelope"); + return; + } + + /* find out how many tracks were actually found... many are zero + \todo is this necessary? */ + for(i = 0, k=0; i < pSmsData->nTracks; i++) + { + if(pSmsData->pFSinFreq[i] > 0.00001) { - if(pSmsData->pFSinFreq[i] > 0.00001) - { - if(pSpecEnvParams->iAnchor != 0) - { - if(k == 0) /* add anchor at beginning */ + if(pSpecEnvParams->iAnchor != 0) + { + if(k == 0) /* add anchor at beginning */ - { - pFreqBuff[k] = 0.0; - pMagBuff[k] = pSmsData->pFSinAmp[i]; - k++; - } - } - pFreqBuff[k] = pSmsData->pFSinFreq[i]; - pMagBuff[k] = pSmsData->pFSinAmp[i]; - k++; + { + pFreqBuff[k] = 0.0; + pMagBuff[k] = pSmsData->pFSinAmp[i]; + k++; } + } + pFreqBuff[k] = pSmsData->pFSinFreq[i]; + pMagBuff[k] = pSmsData->pFSinAmp[i]; + k++; } - /* \todo see if adding an anchor at the max freq helps */ - + } + /* \todo see if adding an anchor at the max freq helps */ - if(k < 1) // how few can this be? try out a few in python - return; - sms_dCepstrum(sizeCepstrum, pSmsData->pSpecEnv, k, pFreqBuff, pMagBuff, - pSpecEnvParams->fLambda, pSpecEnvParams->iMaxFreq); - if(pSpecEnvParams->iType == SMS_ENV_FBINS) - { - sms_dCepstrumEnvelope(sizeCepstrum, pSmsData->pSpecEnv, - pSpecEnvParams->nCoeff, pSmsData->pSpecEnv); - } + if(k < 1) // how few can this be? try out a few in python + return; + sms_dCepstrum(sizeCepstrum, pSmsData->pSpecEnv, k, pFreqBuff, pMagBuff, + pSpecEnvParams->fLambda, pSpecEnvParams->iMaxFreq); + + if(pSpecEnvParams->iType == SMS_ENV_FBINS) + { + sms_dCepstrumEnvelope(sizeCepstrum, pSmsData->pSpecEnv, + pSpecEnvParams->nCoeff, pSmsData->pSpecEnv); + } } |