Реальный выход FFT

Question

Я реализовал fft в ucontroller серии at32ucb, используя библиотеку kiss fft и в настоящее время борюсь с выходом fft. Мое намерение состоит в том, чтобы проанализировать звук, исходящий из пьезо-динамика. В настоящее время частота эхолота составляет 420 Гц, которую я успешно получил с выхода БПФ (перекрестная проверка с помощью осциллографа). Тем не менее, выходная частота составляет лишь половину ожидаемой, если я добавлю сигнал формы функционального генератора в систему. Я подозреваю, что это формула вычисления частотного интервала, которую я ошибся; в настоящее время используется, fft_peak_magnitude_index* частота дискретизации / fft_size. Мой вклад реален и делает реальное БПФ. (output samples = N/2) А также выполняет фильтрацию и управление iir до fft. Любое предложение будет отличной помощью!

            // IIR filter calculation, n = 256 fft points       
        for (ctr=0; ctr<n; ctr++)
        {       
            // filter calculation
            y[ctr] = num_coef[0]*x[ctr];
            y[ctr] += (num_coef[1]*x[ctr-1]) - (den_coef[1]*y[ctr-1]);
            y[ctr] += (num_coef[2]*x[ctr-2]) - (den_coef[2]*y[ctr-2]);
            y1[ctr] = y[ctr] - 510; //eliminate dc offset

            // hamming window
            hamming[ctr] = (0.54-((0.46) * cos(2*M_PI*ctr/n)));
            window[ctr] = hamming[ctr]*y1[ctr];

            fft_input[ctr].r = window[ctr];
            fft_input[ctr].i = 0;
            fft_output[ctr].r = 0;
            fft_output[ctr].i = 0;
        }


        kiss_fftr_cfg fftConfig = kiss_fftr_alloc(n,0,NULL,NULL);
        kiss_fftr(fftConfig, (kiss_fft_scalar * )fft_input, fft_output);

        peak = 0;
        freq_bin = 0;       
        for (ctr=0; ctr<n1; ctr++)
            {   
                fft_mag[ctr] = 10*(sqrt((fft_output[ctr].r * fft_output[ctr].r) + (fft_output[ctr].i * fft_output[ctr].i)))/(0.5*n);

                if(fft_mag[ctr] > peak)
                {
                    peak = fft_mag[ctr];
                    freq_bin = ctr;
                }

            frequency = (freq_bin*(10989/n)); // 10989 is the sampling freq
                //************************************
                //Usart write
                char filtResult[10];
                //sprintf(filtResult, "%04d %04d %04d\n", (int)peak, (int)freq_bin, (int)frequency);
                sprintf(filtResult, "%04d  %04d  %04d\n", (int)x[ctr], (int)fft_mag[ctr], (int)frequency);
                char c;
                char *ptr = &filtResult[0];
                do
                {
                    c = *ptr;
                    ptr++;
                    usart_bw_write_char(&AVR32_USART2, (int)c);
                    // sendByte(c);

                } while (c != '\n');
            }

Answer 1

Основная проблема, скорее всего, заключается в том, как вы заявили fft_input, Исходя из вашего предыдущего вопроса, вы распределяете fft_input как массив kiss_fft_cpx, Функция kiss_fftr с другой стороны ожидаем массив скаляров. Приведя входной массив в kiss_fft_scalar с:

kiss_fftr(fftConfig, (kiss_fft_scalar * )fft_input, fft_output);

KissFFT, по сути, видит массив данных с реальной стоимостью, который содержит нули в каждой второй выборке (то, что вы заполнили как мнимые части). Это фактически версия с повышением частоты дискретизации (хотя и без интерполяции) вашего исходного сигнала, то есть сигнала с удвоенной частотой дискретизации (что не учитывается в вашем freq_bin в frequency преобразование). Чтобы это исправить, я предлагаю вам упаковать ваши данные в kiss_fft_scalar массив:

kiss_fft_scalar fft_input[n];
...
for (ctr=0; ctr<n; ctr++)
{       
    ...
    fft_input[ctr] = window[ctr];
    ...
}
kiss_fftr_cfg fftConfig = kiss_fftr_alloc(n,0,NULL,NULL);
kiss_fftr(fftConfig, fft_input, fft_output);

Также обратите внимание, что при поиске пиковой величины вы, вероятно, интересуетесь только последним наибольшим пиком, а не рабочим максимумом. Таким образом, вы можете ограничить цикл только вычислением пика (используя freq_bin вместо ctr в качестве индекса массива в следующем sprintf заявления при необходимости):

for (ctr=0; ctr<n1; ctr++)
{   
  fft_mag[ctr] = 10*(sqrt((fft_output[ctr].r * fft_output[ctr].r) + (fft_output[ctr].i * fft_output[ctr].i)))/(0.5*n);

  if(fft_mag[ctr] > peak)
  {
    peak = fft_mag[ctr];
    freq_bin = ctr;
  }
} // close the loop here before computing "frequency"

Наконец, при вычислении частоты, связанной с ячейкой с наибольшей величиной, необходимо убедиться, что вычисление выполняется с использованием арифметики с плавающей запятой. Если как я подозреваю n целое число, ваша формула будет выполнять 10989/n фактор с использованием целочисленной арифметики, приводящий к усечению. Это можно просто исправить с помощью:

frequency = (freq_bin*(10989.0/n)); // 10989 is the sampling freq