www.pudn.com > DSP.rar > LLAEDSP.cpp
//--------------------------------------------------------- // Copyright (C) 2005--2006. // All Rights Reserved. //--------------------------------------------------------- // Low level C-like procedure for signal processing #include "LLAEDSP.h" #include#include //--------------------------------------------------------- #define HIFREQ 800 #define LOFREQ 80 //--------------------------------------------------------- void xft(float * vr, float * vi, int order, int s); void shuffle(float * vr, float * vi, int order); int Order_Log2(int num); bool mfcc(float * dp_fft, float * dp_mfcc, float d_sam, int l_nfft, int l_nmfcc); //--------------------------------------------------------- // fft: keeps preemphasized vector of time samples (variable vet_re) // and returns the real (vet_re) and the immaginary (vet_im) parts of its fft //--------------------------------------------------------- void Fft_Call ( float * vet_re, float * vet_im, int vet_dim ) { int one = 1, lx, i, order; float x; float * vr; float * vi; order = Order_Log2(vet_dim); vr = &(vet_re[0]); vr = vr - 1; vi = &(vet_im[0]); vi = vi - 1; xft(vr, vi, order, one); x = (float)(pow(2.0, (float)order)); lx = (int)x; x = (float)sqrt(x); for (i = 1; i <= lx; i++) { *(vr + i) = *(vr + i) / x; *(vi + i) = *(vi + i) / x; } return; } //--------------------------------------------------------- void Mfcc_Call ( float * mfc_vet, int mfc_vet_dim, float * fft_module, int fft_dim, float d_sam ) { float * dp_fft; float * df_mfcc; int l_nfft, l_nmfcc; l_nfft = fft_dim; l_nmfcc = mfc_vet_dim; df_mfcc = &(mfc_vet[0]); df_mfcc = df_mfcc - 1; dp_fft = &(fft_module[0]); dp_fft = dp_fft - 1; mfcc(dp_fft, df_mfcc, d_sam, l_nfft, l_nmfcc); return; } //--------------------------------------------------------- void Pitch_Stress_Call ( float & fCurrFramePitch, float & fCurrFrameStress, float * MagDataP, int iWindowSize, int iSampleRate ) { float * pR, * pI, tmp_Mag; float * pReSpectrumData, * pImSpectrumData; float fSpectrumPeak[2], power, vowelpower, fStress; int iPeakPosition[2]; int i, iSpectrumIndex, iSpectrumSize, iFFTSize; pReSpectrumData = MagDataP; iSpectrumSize = iWindowSize / 2; iFFTSize = iWindowSize; pImSpectrumData = new float [iWindowSize]; power = 0; vowelpower = 0; pR = pReSpectrumData; pR++; for (i = 1; i < iSpectrumSize; i++) { power += (*pR) * (*pR); if (i >= (LOFREQ * iFFTSize / iSampleRate) && i <= (HIFREQ * iFFTSize / iSampleRate)) vowelpower += (*pR) * (*pR); pR++; } memset(pImSpectrumData, 0, iWindowSize * sizeof(float)); Fft_Call(pReSpectrumData, pImSpectrumData, iFFTSize); pR = pReSpectrumData; pI = pImSpectrumData; for (i = 0; i < iFFTSize; i++) { pR[i] = pR[i] * pR[i] + pI[i] * pI[i]; pR[i] = (float)sqrt(pR[i]); } // Extract peak of spectrum in Pitch Window pR = &pReSpectrumData[(iSampleRate / HIFREQ)]; for (i = 0; i < 2; i++) { fSpectrumPeak[i] = pR[0]; iPeakPosition[i] = iSampleRate / HIFREQ; } for (i = 0; i <= (iSampleRate / LOFREQ); i++) { if (pR[1] >= pR[0] && pR[1] >= pR[2]) { if (pR[1] > fSpectrumPeak[0]) { fSpectrumPeak[1] = fSpectrumPeak[0]; iPeakPosition[1] = iPeakPosition[0]; fSpectrumPeak[0] = pR[1]; iPeakPosition[0] = i + iSampleRate / HIFREQ; } else if (pR[1] > fSpectrumPeak[1]) { fSpectrumPeak[1] = pR[1]; iPeakPosition[1] = i + iSampleRate / HIFREQ; } } pR++; } // Stress fStress = vowelpower / ((HIFREQ - LOFREQ) * iFFTSize / iSampleRate); if (fStress < 1) fCurrFrameStress = 0; else fCurrFrameStress = (float)(10 * log10(fStress)); // Pitch tmp_Mag = 1.0; iSpectrumIndex = iFFTSize / iPeakPosition[0]; tmp_Mag = pReSpectrumData[iSpectrumIndex]; fCurrFramePitch = (float)iSpectrumIndex * (float)iSampleRate / (float)iFFTSize; for (i = 0; i < 4; i++) { iSpectrumIndex /= 2; if (iSpectrumIndex < (LOFREQ * iFFTSize / iSampleRate)) break; if (pReSpectrumData[iSpectrumIndex] > tmp_Mag) { tmp_Mag = pReSpectrumData[iSpectrumIndex]; fCurrFramePitch = (float)iSpectrumIndex * (float)iSampleRate / (float)iFFTSize; } } // if ((vowelpower / power) < 0.10) fCurrFramePitch = 0; if (pImSpectrumData) delete [] pImSpectrumData; return; } //--------------------------------------------------------- // ifft: inverse FFT transform. Gets a vector of DFT samples in vet_re // and returns the inverse transform in vet_re and vet_im // Attention: the form of the input array must be specified //--------------------------------------------------------- void Ifft_Call(float * vet_re, float * vet_im, int vet_dim) { int nega = -1, order; float x; int i, lx; float * vr; float * vi; order = Order_Log2(vet_dim); vr = &(vet_re[0]); vr = vr - 1; vi = &(vet_im[0]); vi = vi - 1; x = (float)(pow(2.0, (float)order)); lx = (int)x; x = (float)sqrt((float)x); xft(vr, vi, order, nega); for (i = 1L; i <= lx; i++) { *(vr + i) = *(vr + i) / x; *(vi + i) = *(vi + i) / x; } return; } //--------------------------------------------------------- int Order_Log2(int num) { int act_num, order; order = 0; do { order++; act_num = 1 << order; if (act_num > num) return -act_num; } while (act_num != num); return order; } //--------------------------------------------------------- // given the complex vector with real part xr[1..n] and imaginary part, // xi[1..n], computes the FFT on NN=2^order points and puts results in the // xr[1..NN], xi[1..NN] vector position, if SIGN>=0 // computes the reverse FT if SIGN is negative // the routine SHUFFLE is referenced. //--------------------------------------------------------- void xft(float * vr, float * vi, int order, int s) { register int i, j, n, n2; int z, k, i_or; float wr, wi, arg, dwr, dwi, tr, ti, theta; if (s >= 0) theta = (float)acos(-1.0); else theta = (float)(-1.0 * acos(-1.0)); shuffle(vr, vi, order); n2 = 1; i_or = (int)order; for (i = 1; i <= i_or; i++) { n = n2; n2 += n2; wr = (float)1.0; wi = (float)0.0; arg = theta / (float)n; dwr = (float)cos((float)arg); dwi = (float)sin((float)arg); for (j = 1; j <= n; j++) { for (k = j; k <= (1 << i_or); k += n2) { z = k + n; tr = vr[z] * wr + vi[z] * wi; ti = vi[z] * wr - vr[z] * wi; vr[z] = vr[k] - tr; vi[z] = vi[k] - ti; vr[k] += tr; vi[k] += ti; } tr = wr * dwr - wi * dwi; wi = wi * dwr + wr * dwi; wr = tr; } } return; } //--------------------------------------------------------- void shuffle(float * vr, float * vi, int order) { int i, j, k, n2, i_or; float t; i_or = (int)order; n2 = (1 << i_or) >> 1; j = 1; for (i = 1; i < (1 << i_or); i++) { if (i < j) { t = vr[i]; vr[i] = vr[j]; vr[j] = t; t = vi[i]; vi[i] = vi[j]; vi[j] = t; } k = n2; while (k < j) { j -= k; k >>= 1; } j += k; } return; } //--------------------------------------------------------- // REFERENCES // Davis, Mermelstein "Comparison of parametric representation ..." // ASSP Vol.28, No.4, August 80, pp.357 //--------------------------------------------------------- bool mfcc(float * dp_fft, float * dp_mfcc, float d_sam, int l_nfft, int l_nmfcc) { char * buffer; float band, rowamp, onek; float * filterbank; int i, j, n, k, nfilterbank; int delta, left, center, rigth; float a, b, pigr; float zero = 0.0; float fdelta; float factor = (float)1.2; band = d_sam / (float)2.0; rowamp = band / (float)l_nfft; if (l_nfft != 64 && l_nfft != 128 && l_nfft != 256 && l_nfft != 512 && l_nfft != 1024) return false; buffer = new char [sizeof(float) * (l_nfft + 4)]; if (!buffer) return false; filterbank = (float *)&(buffer[0]); for (i = 0; i <= l_nfft; i++) filterbank[i] = 0; //******************************************************** //* * //* INITIALIZE FILTER DEFINITION IN LINEAR PART * //* * //******************************************************** i = 1; onek = (float)1000.0 / rowamp; // rows in 1kHz fdelta = onek / (float)10.0; // we wants about 10 cep uder 1kHz delta = (int)fdelta; n = 1; left = n; center = n + delta; rigth = n + 2 * delta; //******************************************************** //* * //* NOW COMPUTES THE LINEAR SCALE BANKFILTER f<1kHz * //* * //******************************************************** // until the central boundary is lower than 1kHz while (n < (int)onek) { filterbank[i] = 0.0; // compute the decrement step unit for triangular filters b = (float)2.0 + (float)rigth - (float)left; b = (b - (float)1.0) / (float)2.0; b = (float)1.0 / b; // b is the step unit for (j = left; j <= rigth; j++ ) { // a is the point distance from the centre of the // triangular filter a = (float)fabs((float)center - (float)j); a = (float)1.0 - a * b; // a is the factor<1 to amplify row(j) // update filter i in filterbank computation filterbank[i] = filterbank[i] + dp_fft[j] * a; } // NOW // 1) compute next filter i -> i+1 // 2) delta is still the same (linear grows) // 3) filter center moves up of delta points // 4) left and rigth boundary are updated i = i + 1; n = n + delta; left = n; center = n + delta; rigth = n + 2 * delta; } //******************************************************** //* * //* NOW COMPUTES THE LOG SCALE BANKFILTER f>1kHz * //* * //******************************************************** // until the rigth boundary is inside the FFT while (rigth < l_nfft) { filterbank[i] = 0.0; // compute the decrement step unit for triangular filters b = (float)2.0 + (float)rigth - (float)left; b = (b - (float)1.0) / (float)2.0; b = (float)1.0 / b; // b is the step unit for (j = left; j <= rigth; j++) { // a is the point distance from the centre of the // triangular filter a = (float)fabs((float)center - (float)j); a = (float)1.0 - a * b; // a is the factor<1 to amplify row(j) // update filter i in filterbank computation filterbank[i] = filterbank[i] + dp_fft[j] * a; } // NOW // 1) compute next filter i -> i+1 // 2) delta grows exponentially delta -> delta*factor // 3) filter center moves up of delta points // 4) left and rigth boundary are updated i = i +1; fdelta = factor * fdelta; delta = (int) fdelta; n = n + delta; left = n; center = n + delta; rigth = n + 2 * delta; } //******************************************************** //* * //* LOG ENERGY FILTERBANK AND CEPSTRUM COMPUTATION * //* * //******************************************************** nfilterbank = i - 1; // the number of filter if (nfilterbank < l_nmfcc) return false; // LOG computation of the filterbank, // no negative values are allowed, if E<1 => E=1. for (i = 1; i <= nfilterbank; i++) { if (filterbank[i] >= 1.0) { filterbank[i] = (float)(20.0 * log10(filterbank[i])); } else { filterbank[i] = 0.0; } } // the cosine transformation to obtain MFCC pigr = (float)acos(-1.0); for (i = 1; i <= l_nmfcc; i++) { dp_mfcc[i] = (float)0.0; for (k = 1; k <= nfilterbank; k++) { a = (float)k - (float)0.5; a = a * pigr * (float)i; a = a / (float)nfilterbank; a = (float)cos(a); dp_mfcc[i] += filterbank[k] * a; } } if (buffer) delete buffer; return true; }