MFCC

пожалуйста, дайте формулы по которым считается MFCC.

формула обратного косинусоидального преобразования тут
расскажу в двух словах один из способов, как они высчитываются
1) Сигнал переводится в спектральную область с помощью БПФ (быстрого преобразоавания Фурье)
2) Спектральные коэффициенты усредняются с использованием параметра количества каналов фильтра, которые представляют из себя гребенку треугольных окон, расположенных на мел-шкале (обычно треугольные фильтры перекрывают область частот от нуля и вплоть до частоты Найквиста ).
3) Мел-частотные кепстральные коэффициенты вычисляются при помощи обратного косинусного преобразования (см. формулу по ссылке).

Сообщение отредактировано: goodass - 29.04.08, 10:15

спасибо

Хочу протестировать правильность нахождения мел-частотных кепстральных коэффициентов в моей программе.

Может кто-нибудь выложить кусок звукового сигнала с выделенными из него MFCC?

Я Вам настойчиво советую внимательно прочитать предыдущую ветку на эту тему:

MFCC how????

я не смог разобраться в логических переходах между сообщениями в этой ветке.

если я правильно понял, http://htk.eng.cam.ac.uk/ - на этом сайте я могу скачать звуковые фрагменты с выделенным MFCC?

Хм, ну попробую донести мысль всё-таки.

Вычисленные коэффициенты MFCC зависят от большого числа параметров - числа фильтров, формы фильтров, числа точек в спектре, основания логарифма и т.д. Поэтому с вероятностью 100% две различные реализации MFCC дадут совершенно различные числа. Поэтому рекомендуется перейти по ссылке http://htk.eng.cam.ac.uk, зарегистрироваться там, и скачать уже скомпилированные под Windows приложения, в частности вычисляющие MFCC. Документация к HTK на русском языке доступна здесь:

http://aprodeus.narod.ru/mmedia.htm#Rech

Используя HCopy и HList для вычисления коэффициентов и сверяясь с исходым кодом HTK возможно проверить работу приложения.

спасибо огромное

Доброго всем дня.
Вопрос знатокам HTK по MFCC:
код написанный с использованием IPP и HList выдают разные результаты, с чем это может быть связано?
IPP:
01: -29.003057
02: -20.139099
03: -4.3270502
04: -13.346766
05: -18.708612
06: -5.8949571
07: 31.116928
08: 27.441202
09: 21.219500
10: 39.831188
11: 13.616776
12: -3.3898118
E: 9.6506577

HList:
------------------------------------------------- Source: ru_0001.wav --------------------------------------------------
Sample Bytes: 2 Sample Kind: WAVEFORM
Num Comps: 1 Sample Period: 62.5 us
Num Samples: 257278 File Format: WAV
-------------------------------------------------------- Target --------------------------------------------------------
Sample Bytes: 52 Sample Kind: MFCC_E
Num Comps: 13 Sample Period: 10000.0 us
Num Samples: 1606 File Format: HTK
------------------------------------------------ Observation Structure -------------------------------------------------
x: MFCC-1 MFCC-2 MFCC-3 MFCC-4 MFCC-5 MFCC-6 MFCC-7 MFCC-8 MFCC-9 MFCC-10 MFCC-11 MFCC-12 E
---------------------------------------------------- Samples: 0->0 -----------------------------------------------------
0: -10.925 -10.630 -1.796 -8.355 -11.068 -5.632 19.115 14.508 12.260 21.412 7.453 -1.746 9.651
--------------------------------------------------------- END ----------------------------------------------------------

config:
# Coding parameters
SOURCEFORMAT = WAV
TARGETKIND = MFCC_E
TARGETRATE = 100000.0
WINDOWSIZE = 250000.0
USEHAMMING = T
PREEMCOEF = 0.97
NUMCHANS = 24
CEPLIFTER = 22
NUMCEPS = 12
ENORMALISE = F
LOFREQ = 100
HIFREQ = 4000
ESCALE = 1.0
SILFLOOR = 50.0
USEPOWER = T

код:

    /* Initialize the structures */
    ippsMelFBankInitAlloc_32f(&fbank, /* return the structure pointer */
        &fft_order, /* return the FFT length */
        400,    /* 25ms window/512 point FFT */
        16000,  /* sample rate */
        100,    /* lowest frequency of interest */
        4000,   /* highest frequency of interest */
        24,     /* number of filter banks */
        1127, /* mel-scale factor 1 */
        700,    /* mel-scale factor 2 */
        (IppMelMode)(IPP_FBANK_MELWGT | IPP_POWER_SPECTRUM));
    ippsDCTLifterInitAlloc_32f(&dctl, /* return the structure pointer */
        24, /* filter bank channels */
        12, /* number of MFCC coefficients */
        22, /* liftering */
        1.0f); /* no scaling */
    fft_len=1<<fft_order;
    frame_buffer=ippsMalloc_32f(fft_len);
    fbank_buffer=ippsMalloc_32f(24);
    tempmfccs=ippsMalloc_32f(13*(sample_num-400+160)/160);
    mfcc_cur=tempmfccs;
    /* Calculate MFCC features */
    for (i=j=0; i+WindowSize<=sample_num; i+=StepSize,mfcc_cur+=13,j++) {
        /* Organize the input wave data into a frame */
        ippsCopyWithPadding_32f(&samples[i],400,frame_buffer,fft_len);
        ippsDotProd_32f(frame_buffer,frame_buffer,fft_len,&LogEnergy);
        /* Pre-emphasize the input signal with factor 0.97 */
        ippsPreemphasize_32f(frame_buffer,400,0.97f);
        frame_buffer[0]*=(1.0f-0.97f);
        /* Add the hamming window to the input signal */
        ippsWinHamming_32f_I(frame_buffer,WindowSize);
        /* Perform the filter bank analysis */
        ippsEvalFBank_32f(frame_buffer,fbank_buffer,fbank);
        ippsThreshold_LTVal_32f_I(fbank_buffer, 24, 1.0f, 1.0f);
        ippsLn_32f_I(fbank_buffer, 24);
        /* Perform the DCT analysis and liftering */
        ippsDCTLifter_32f(fbank_buffer,mfcc_cur,dctl);
        mfcc_cur[12] = (float) log ((double) LogEnergy);
    }
    /* Normalize log energy */
    //ippsNormEnergy_32f(tempmfccs+12,13,(sample_num-400+160)/160,50.0f,1.0f);

Нашёл на просторах интернета вот такой класс MFCC!
Каково мнение гуру?
По идее,я посмотрел,там есть всё!И БПФ,и окно Хемминга,и логарифмирование и DCT и на выходе массив кепстров!На входе должен быть массив double[].Я пока провожу "эксперимент"...

using System;
using System.Collections.Generic;
using System.Text;
 
namespace speakerRecognition
{
    class Mfcc
    {
        static int my_rint(double x)
        {
            if (2 * x == (double)Math.Round(2 * x))
                x += 0.0001;
 
            return ((int)Math.Round(x));
        }
 
        static public double[] mfcc(ref double[] signal)
        {
            double lowestFrequency = 133.3333;
            double linearFilters = 13;
            double linearSpacing = 66.66666666;
            double logFilters = 27;
            double logSpacing = 1.0711703;
            double fftSize = 512;
            double cepstralCoefficients = 13;
            double windowSize = 256;
            double samplingRate = 44100;
            double totalFilters;
            double[] mfccFilterWeights;
            double[] freqs;
            double[] triangleHeight;
            double[] lower;
            double[] upper;
            double[] center;
            double[] fftFreqs;
            double[] hamWindow;
            double[] mfccDCTMatrix;
            double[] preEmphasized;
            double first;
            double last;
            double[] fftData;
            double[] empreinte;
            double maxi;
            double[] earMag;
            double[] ceps;
            double loga;
            double windowStep;
            double cols;
            double frameRate = 100;
            double[] result;
 
            int i, j, k;
 
 
            preEmphasized = new double[1024];
 
            totalFilters = linearFilters + logFilters;
            freqs = new double[(int)totalFilters + 2];
            for (i = 0; i < linearFilters; i++)
                freqs[i] = lowestFrequency + i * linearSpacing;
            for (i = (int)linearFilters; i < (int)totalFilters + 2; i++)
                freqs[i] = freqs[(int)linearFilters - 1] * Math.Pow(logSpacing, i - linearFilters + 1);
            lower = new double[(int)totalFilters];
            upper = new double[(int)totalFilters];
            center = new double[(int)totalFilters];
            triangleHeight = new double[(int)totalFilters];
            fftData = new double[(int)fftSize];
 
            for (i = 0; i < (int)totalFilters; i++)
                lower[i] = freqs[i];
            for (i = 1; i < (int)totalFilters + 1; i++)
                center[i - 1] = freqs[i];
            for (i = 2; i < (int)totalFilters + 2; i++)
                upper[i - 2] = freqs[i];
            hamWindow = new double[(int)windowSize];
            mfccFilterWeights = new double [(int)(totalFilters * fftSize)];
            for (i = 0; i < totalFilters; i++)
                triangleHeight[i] = 2 / (upper[i] - lower[i]);
            fftFreqs = new double[(int)fftSize];
            for (i = 0; i < fftSize; i++)
                fftFreqs[i] = (i / fftSize) * samplingRate;
 
            for (i = 0; i < totalFilters; i++)
            {
                for (j = 0; j < (int)fftSize; j++)
                {
                    if ((fftFreqs[j] > lower[i]) && (fftFreqs[j] <= center[i]))
                        mfccFilterWeights[(int)fftSize * i + j] = triangleHeight[i] * (fftFreqs[j] - lower[i]) / (center[i] - lower[i]);
                    if ((fftFreqs[j] > center[i]) && (fftFreqs[j] < upper[i]))
                        mfccFilterWeights[(int)fftSize * i + j] = (triangleHeight[i] * (fftFreqs[j] - lower[i]) / (center[i] - lower[i])) + (triangleHeight[i] * (upper[i] - fftFreqs[j]) / (upper[i] - center[i]));
                }
            }
 
            for (i = 0; i < windowSize; i++)
                hamWindow[i] = 0.54 - 0.46 * Math.Cos(2 * Math.PI * i / windowSize);
            mfccDCTMatrix = new double[(int)(cepstralCoefficients * totalFilters)];
            for (i = 0; i < cepstralCoefficients; i++)
                for (j = 0; j < totalFilters; j++)
                    mfccDCTMatrix[j * (int)cepstralCoefficients + i] = 1 / Math.Sqrt(totalFilters / 2) * Math.Cos(i * (2 * j + 1) * Math.PI / 2 / totalFilters);
            for (i = 0; i < totalFilters; i++)
                mfccDCTMatrix[i * (int)cepstralCoefficients] *= 1 / Math.Sqrt(2);
            for (i = 1; i < 1024; i++)
                preEmphasized[i] = signal[i] + signal[i - 1] * -0.97;
            preEmphasized[0] = signal[0];
 
 
            earMag = new double[(int)totalFilters];
            windowStep = samplingRate / frameRate;
            cols = (int)((1024 - windowSize) / windowStep);
            ceps = new double[(int)(sizeof(double) * cols * linearFilters)];
            for (i = 0; i < cols; i++)
            {
                first = i * windowStep + 1;
                last = first + windowSize - 1;
                for (j = 0; j < fftSize; j++)
                    fftData[j] = 0;
                for (j = 0; j < windowSize; j++)
                    fftData[j] = preEmphasized[(int)my_rint(first + j - 1)] *  hamWindow[j];
                double[] fftMag = FourierTransform.FFT(ref signal);
                for (j = 0; j < totalFilters; j++)
                {
                    loga = 0;
                    for (k = 0; k < fftSize; k++)
                        loga += fftMag[k] * mfccFilterWeights[k + j * (int)fftSize];
 
                    earMag[j] = Math.Log10(loga);
                }
                for (j = 0; j < linearFilters; j++)
                {
                    loga = 0;
                    for (k = 0; k < totalFilters; k++)
                        loga +=
                          earMag[k] * mfccDCTMatrix[j + k * (int)linearFilters];
                    ceps[j + i * (int)linearFilters] = loga;
                }
            }
 
 
            empreinte = new double[12];
 
            for (i = 1; i < 13; i++)
                empreinte[i - 1] = (ceps[i] + ceps[i + 13] + ceps[i + 26] + ceps[i + 39] + ceps[i + 52] + ceps[i + 65]);
 
            maxi = 0;
            for (i = 0; i < 12; i++)
            {
                if (maxi * maxi < empreinte[i] * empreinte[i])
                    maxi = empreinte[i];  
            }
            for (i = 0; i < 12; i++) empreinte[i] = (empreinte[i] / Math.Abs(maxi)) * 2;
            
            return empreinte;
 
        }
    }
}

Сообщение отредактировано: Dj_SheLL - 10.04.11, 12:14

Насчёт БПФ.
Там вроде и прямое и обратное можно.
Брал вот этот модуль:

/// <license>
/// This is a port of the SciMark2a Java Benchmark to C# by
/// Chris Re (cmr28@cornell.edu) and Werner Vogels (vogels@cs.cornell.edu)
///
/// For details on the original authors see http://math.nist.gov/scimark2
///
/// This software is likely to burn your processor, bitflip your memory chips
/// anihilate your screen and corrupt all your disks, so you it at your
/// own risk.
/// </license>
 
using System;
 
namespace SciMark2
{
    
    /// <summary>Computes FFT's of complex, double precision data where n is an integer power of 2.
    /// This appears to be slower than the Radix2 method,
    /// but the code is smaller and simpler, and it requires no extra storage.
    /// </P>
    /// </summary>
    ///
    /// <author>
    /// Bruce R. Miller bruce.miller@nist.gov,
    /// Derived from GSL (Gnu Scientific Library),
    /// GSL's FFT Code by Brian Gough bjg@vvv.lanl.gov
    /// </author>
    
    public class FFT
    {
        
        public static double num_flops(int N)
        {
            double Nd = (double) N;
            double logN = (double) log2(N);
            
            return (5.0 * Nd - 2) * logN + 2 * (Nd + 1);
        }
        
        
        /// <summary>
        /// Compute Fast Fourier Transform of (complex) data, in place.
        /// </summary>
        public static void  transform(double[] data)
        {
            transform_internal(data, - 1);
        }
        
        /// <summary>
        /// Compute Inverse Fast Fourier Transform of (complex) data, in place.
        /// </summary>
        public static void  inverse(double[] data)
        {
            transform_internal(data, + 1);
            // Normalize
            int nd = data.Length;
            int n = nd / 2;
            double norm = 1 / ((double) n);
             for (int i = 0; i < nd; i++)
                data[i] *= norm;
        }
        
        /// <summary>
        /// Accuracy check on FFT of data. Make a copy of data, Compute the FFT, then
        /// the inverse and compare to the original.  Returns the rms difference.
        /// </summary>
        public static double test(double[] data)
        {
            int nd = data.Length;
            // Make duplicate for comparison
            double[] copy = new double[nd];
            Array.Copy(data, 0, copy, 0, nd);
            // Transform & invert
            transform(data);
            inverse(data);
            // Compute RMS difference.
            double diff = 0.0;
             for (int i = 0; i < nd; i++)
            {
                double d = data[i] - copy[i];
                diff += d * d;
            }
            return Math.Sqrt(diff / nd);
        }
        
        /// <summary>
        /// Make a random array of n (complex) elements.
        /// </summary>
        public static double[] makeRandom(int n)
        {
            int nd = 2 * n;
            double[] data = new double[nd];
            System.Random r = new System.Random();
            for (int i = 0; i < nd; i++)
                data[i] = r.NextDouble();
            return data;
        }
        
        protected internal static int log2(int n)
        {
            int log = 0;
             for (int k = 1; k < n; k *= 2, log++)
                ;
            if (n != (1 << log))
                throw new ApplicationException("FFT: Data length is not a power of 2!: " + n);
            return log;
        }
        
        protected internal static void  transform_internal(double[] data, int direction)
        {
            if (data.Length == 0)
                return ;
            int n = data.Length / 2;
            if (n == 1)
                return ;
            // Identity operation!
            int logn = log2(n);
            
            /* bit reverse the input data for decimation in time algorithm */
            bitreverse(data);
            
            /* apply fft recursion */
            /* this loop executed log2(N) times */
             for (int bit = 0, dual = 1; bit < logn; bit++, dual *= 2)
            {
                double w_real = 1.0;
                double w_imag = 0.0;
                
                double theta = 2.0 * direction * Math.PI / (2.0 * (double) dual);
                double s = Math.Sin(theta);
                double t = Math.Sin(theta / 2.0);
                double s2 = 2.0 * t * t;
                
                /* a = 0 */
                 for (int b = 0; b < n; b += 2 * dual)
                {
                    int i = 2 * b;
                    int j = 2 * (b + dual);
                    
                    double wd_real = data[j];
                    double wd_imag = data[j + 1];
                    
                    data[j] = data[i] - wd_real;
                    data[j + 1] = data[i + 1] - wd_imag;
                    data[i] += wd_real;
                    data[i + 1] += wd_imag;
                }
                
                /* a = 1 .. (dual-1) */
                 for (int a = 1; a < dual; a++)
                {
                    /* trignometric recurrence for w-> exp(i theta) w */
                    {
                        double tmp_real = w_real - s * w_imag - s2 * w_real;
                        double tmp_imag = w_imag + s * w_real - s2 * w_imag;
                        w_real = tmp_real;
                        w_imag = tmp_imag;
                    }
                     for (int b = 0; b < n; b += 2 * dual)
                    {
                        int i = 2 * (b + a);
                        int j = 2 * (b + a + dual);
                        
                        double z1_real = data[j];
                        double z1_imag = data[j + 1];
                        
                        double wd_real = w_real * z1_real - w_imag * z1_imag;
                        double wd_imag = w_real * z1_imag + w_imag * z1_real;
                        
                        data[j] = data[i] - wd_real;
                        data[j + 1] = data[i + 1] - wd_imag;
                        data[i] += wd_real;
                        data[i + 1] += wd_imag;
                    }
                }
            }
        }
        
        
        protected internal static void  bitreverse(double[] data)
        {
            /* This is the Goldrader bit-reversal algorithm */
            int n = data.Length / 2;
            int nm1 = n - 1;
            int i = 0;
            int j = 0;
             for (; i < nm1; i++)
            {
                
                //int ii = 2*i;
                int ii = i << 1;
                
                //int jj = 2*j;
                int jj = j << 1;
                
                //int k = n / 2 ;
                int k = n >> 1;
                
                if (i < j)
                {
                    double tmp_real = data[ii];
                    double tmp_imag = data[ii + 1];
                    data[ii] = data[jj];
                    data[ii + 1] = data[jj + 1];
                    data[jj] = tmp_real;
                    data[jj + 1] = tmp_imag;
                }
                
                while (k <= j)
                {
                    //j = j - k ;
                    j -= k;
                    
                    //k = k / 2 ;
                    k >>= 1;
                }
                j += k;
            }
        }
    }
}

Сообщение отредактировано: Dj_SheLL - 11.04.11, 07:27

Dj_SheLL, спасибо за информацию. Поделись пожалуйста впечатлениями от использования класса MFCC. Просто я начал писать свой класс, но столкнулся с некоторыми сложностями, так что возможно буду использовать тот, который ты нашел. Как вообще проверить, правильные ли результаты выдает MFCC или нет?

Внешне всё вроде реализованно верно,НО!
Я так и не добился получения кепстров-у меня при вычислениях +/- "бесконечность" получается.Проблемы с длинной чисел.
Если что, данные я получал с буфера с помощью NAudio.По идее с других библиотек данные не должны отличаться ни чем!
Так что надо допиливать класс под "народные нужды",чтобы всем пригодилось...либо пользоваться другими инструментами!

Насчет +/- "бесконечности" ничего не могу посоветовать) Я проверил алгоритм и сравнил с теми, которые реализованы на других языках и лежат на просторах интернета. В целом он не вызывает подозрений.
Что касается кепстров, то у меня они получились в интервале [2;-1]. NAudio я не использовал.
У меня вопрос: перед тем, как записывать звук, мне нужно указать AudioFrameSize. Я так понимаю, что в этот параметр записывается сколько милисекунд будет длится каждый сэмпл. Вляет ли этот параметр на значение кепстров и если да, то сколько мс там указывать?

Цитата Lohoboishe @ 14.04.11, 11:57

Что касается кепстров, то у меня они получились в интервале [2;-1]

А вообще какой диапозон должен быть?

Цитата Lohoboishe @ 14.04.11, 11:57

У меня вопрос: перед тем, как записывать звук, мне нужно указать AudioFrameSize. Я так понимаю, что в этот параметр записывается сколько милисекунд будет длится каждый сэмпл. Вляет ли этот параметр на значение кепстров и если да, то сколько мс там указывать?

Ну все берут от 20 до 30.Насчёт влияния-не знаю.