Что может привести к тому, что данные БПФ будут иметь пики на неправильных частотах?

Question

Я реализую функцию определения высоты тона FFT на iPhone, используя Apple Accelerate Framework, как обсуждалось много раз здесь раньше.

Я понимаю сдвиги фаз, частоты дискретизации и исследовал несколько тюнеров с открытым исходным кодом, которые используют методы FFT (простое определение высоты тона, автокорреляция, кепстр и т. П.) Для определения высоты тона. Вот моя проблема:

Мои результаты БПФ постоянно отключаются на 5-10 Гц (+/-), даже когда бины разнесены всего на 1-2 Гц. Я пробовал разные алгоритмы и даже простое БПФ Выборка с высоким разрешением показывает пики величины в, казалось бы, неправильных местах. Это не последовательное смещение; некоторые слишком высоки, некоторые слишком низки.

Например, тон 440 Гц равен 445,2 Гц; 220 Гц как 214 Гц; 880 Гц как 874 Гц; 1174 Гц как 1183 Гц с использованием тонального генератора. Подобный тюнер с открытым исходным кодом для Mac, использующий почти точно такие же алгоритмы, не имеет проблем с точным определением высоты тона. (Эти различия отличаются, когда на устройстве, чем симулятор, но они по-прежнему выключены.)

Я не думаю, что проблема заключается в разрешении бина, потому что между фактическим тоном и обнаруженным всплеском амплитуды часто бывает несколько бинов. Это как если бы вход просто слышал неправильную высоту.

Я вставил свой код ниже. Общий поток прост:

Нажмите шаг в буфер FFT -> Окно Hann -> FFT -> Phase / Magnitude -> Max pitch is not.

enum {
    kOversample = 4,
    kSamples = MAX_FRAME_LENGTH,
    kSamples2 = kSamples / 2,
    kRange = kSamples * 5 / 16,
    kStep = kSamples / kOversample
};



const int PENDING_LEN = kSamples * 5;
static float pendingAudio[PENDING_LEN * sizeof(float)];
static int pendingAudioLength = 0;

- (void)processBuffer {
    static float window[kSamples];
    static float phase[kRange];
    static float lastPhase[kRange];
    static float phaseDeltas[kRange];
    static float frequencies[kRange];
    static float slidingFFTBuffer[kSamples];
    static float buffer[kSamples];

    static BOOL initialized = NO;
    if (!initialized) {
        memset(lastPhase, 0, kRange * sizeof(float));

        vDSP_hann_window(window, kSamples, 0);
        initialized = YES;
    }

    BOOL canProcessNewStep = YES;
    while (canProcessNewStep) {        

        @synchronized (self) {
            if (pendingAudioLength < kStep) {
                break; // not enough data
            }            
            // Rotate one step's worth of pendingAudio onto the end of slidingFFTBuffer
            memmove(slidingFFTBuffer, slidingFFTBuffer + kStep, (kSamples - kStep) * sizeof(float));
            memmove(slidingFFTBuffer + (kSamples - kStep), pendingAudio, kStep * sizeof(float));
            memmove(pendingAudio, pendingAudio + kStep, (PENDING_LEN - kStep) * sizeof(float));
            pendingAudioLength -= kStep;   
            canProcessNewStep = (pendingAudioLength >= kStep);
        }

        // Hann Windowing
        vDSP_vmul(slidingFFTBuffer, 1, window, 1, buffer, 1, kSamples);      
        vDSP_ctoz((COMPLEX *)buffer, 2, &splitComplex, 1, kSamples2);        

        // Carry out a Forward FFT transform.
        vDSP_fft_zrip(fftSetup, &splitComplex, 1, log2f(kSamples), FFT_FORWARD);        

        // magnitude to decibels
        static float magnitudes[kRange];        
        vDSP_zvmags(&splitComplex, 1, magnitudes, 1, kRange);        
        float zero = 1.0;
        vDSP_vdbcon(magnitudes, 1, &zero, magnitudes, 1, kRange, 0); // to decibels

        // phase
        vDSP_zvphas(&splitComplex, 1, phase, 1, kRange); // compute magnitude and phase        
        vDSP_vsub(lastPhase, 1, phase, 1, phaseDeltas, 1, kRange); // compute phase difference
        memcpy(lastPhase, phase, kRange * sizeof(float)); // save old phase

        double freqPerBin = sampleRate / (double)kSamples;
        double phaseStep = 2.0 * M_PI * (float)kStep / (float)kSamples;

        // process phase difference ( via https://stackoverflow.com/questions/4633203 )
        for (int k = 1; k < kRange; k++) {
            double delta = phaseDeltas[k];
            delta -= k * phaseStep;  // subtract expected phase difference
            delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
            delta /= phaseStep;  // calculate diff from bin center frequency
            frequencies[k] = (k + delta) * freqPerBin;  // calculate the true frequency
        }               

        NSAutoreleasePool *pool = [[NSAutoreleasePool alloc] init];

        MCTunerData *tunerData = [[[MCTunerData alloc] initWithSize:MAX_FRAME_LENGTH] autorelease];        

        double maxMag = -INFINITY;
        float maxFreq = 0;
        for (int i=0; i < kRange; i++) {
            [tunerData addFrequency:frequencies[i] withMagnitude:magnitudes[i]];
            if (magnitudes[i] > maxMag) {
                maxFreq = frequencies[i];
                maxMag = magnitudes[i];
            }
        }

        NSLog(@"Max Frequency: %.1f", maxFreq);

        [tunerData calculate];

        // Update the UI with our newly acquired frequency value.
        [self.delegate frequencyChangedWithValue:[tunerData mainFrequency] data:tunerData];

        [pool drain];
    }

}

OSStatus renderCallback(void *inRefCon, AudioUnitRenderActionFlags *ioActionFlags, 
                       const AudioTimeStamp *inTimeStamp, UInt32 inBusNumber, UInt32 inNumberFrames, 
                       AudioBufferList *ioData)
{
    MCTuner* tuner = (MCTuner *)inRefCon;    

    OSStatus err = AudioUnitRender(tuner->audioUnit, ioActionFlags, inTimeStamp, 1, inNumberFrames, tuner->bufferList);
    if (err < 0) {
        return err;
    }

    // convert SInt16 to float because iOS doesn't support recording floats directly
    SInt16 *inputInts = (SInt16 *)tuner->bufferList->mBuffers[0].mData;

    @synchronized (tuner) {
        if (pendingAudioLength + inNumberFrames < PENDING_LEN) {

            // Append the audio that just came in into the pending audio buffer, converting to float
            // because iOS doesn't support recording floats directly
            for(int i = 0; i < inNumberFrames; i++) {
                pendingAudio[pendingAudioLength + i] = (inputInts[i] + 0.5) / 32767.5;
            }
            pendingAudioLength += inNumberFrames;
        } else {
            // the buffer got too far behind. Don't give any more audio data.
            NSLog(@"Dropping frames...");
        }
        if (pendingAudioLength >= kStep) {
            [tuner performSelectorOnMainThread:@selector(processBuffer) withObject:nil waitUntilDone:NO];
        }
    }

    return noErr;
}

Answer 1

Я не детально рассмотрел ваш код, но это выскочило прямо на меня:

vDSP_zvmags(&splitComplex, 1, magnitudes, 1, kRange);

Важно помнить, что результат от реального к сложному fft упакован внесколько странный макет.Если действительная и мнимая части j-го коэффициента Фурье обозначены через R (j) и I (j), компоненты real и imag объекта splitComplex имеют следующее содержимое:

.real = {  R(0) , R(1), R(2), ... , R(n/2 - 1) } 
.imag = { R(n/2), I(1), I(2), ... , I(n/2 - 1) }

Таким образом, ваш расчет величины делает что-то немного странное;первая запись в вашем векторе величины sqrt(R(0)^2 + R(n/2)^2), где она должна быть |R(0)|.Я не тщательно изучил все константы, но кажется вероятным, что это приводит к частичной ошибке, когда вы теряете полосу Найквиста (R(n/2)) или аналогичную.Подобная ошибка может привести к тому, что полосы частот будут рассматриваться как немного более широкие или более узкие, чем они есть на самом деле, что приведет к небольшому увеличению или уменьшению высоты тона во всем диапазоне, что соответствуетчто ты видишь.

Answer 2

БПФ - это бензопила, а не скальпель. В общем, для проверки FFT-кодирования с проверкой реальности: (1) тестируйте с использованием теоремы Парсеваля (средняя квадратичная амплитуда во временной области должна в пределах округления равна сумме вашего спектра) и (2) обратного БПФ и просто слушайте его. Извините, но вы, кажется, ожидаете слишком много абсолютной точности от БПФ. Вы просто не получите это. Есть, однако, контрольный список мелочей, чтобы проверить в вашем коде. Большинство алгоритмов перемещают DC и Nyquist, чтобы распределить память даже не приходилось, но вы должны вручную переместить термин Nyquist туда, где он принадлежит, и истинно обнулить различные вещи:

A.realp[NOVER2] = A.imagp[0];   // move real Nyquist term to where it belongs
A.imagp[NOVER2] = 0.0;          // this is zero
A.imagp[0] = 0.0;               // make this a true zero

В аудиоданных DC должно быть равно нулю (например, амплитуды имеют среднее значение нуля), но в маленьких окнах это может не быть. Я оставляю это в покое. Вы делаете гораздо больше, чем нужно, чтобы найти максимальную корзину (комментарий о фазовом вокодере правильный). ИМХО использование окна Хамм снижает точность. У меня есть намного лучшие результаты, дополняющие конец реальных данных лотами (4x) нулей. Удачи.

Answer 3

Я уверен, что это не было ничего в моем алгоритме;скорее, что-то не так с моим использованием AUGraph от Apple.Когда я удалил это, чтобы использовать обычный аудиоустройство без настройки графика, я смог заставить его правильно распознавать высоту тона.

Answer 4

Похоже, вы используете не только БПФ, но и фазовый вокодер после БПФ для регулировки предполагаемых частот бина.В зависимости от усиления и пределов фазы, регулировки частоты фазового вокодера могут вывести оценочную частоту далеко за пределы ширины ячейки БПФ.Если это происходит, использование более узких лотков (более длинное БПФ) не поможет.Возможно, вы захотите проверить работоспособность, чтобы увидеть, хотите ли вы, чтобы пиковые частоты выходили за пределы их частотных интервалов БПФ.Или попробуйте убрать фазовый вокодер, чтобы увидеть, возвращает ли одно БПФ более вменяемые результаты.