Делать БПФ в режиме реального времени

Question

Я хочу сделать БПФ аудиосигнала в реальном времени, то есть, когда человек говорит в микрофон.Я получу данные (я делаю это с помощью portaudio, если было бы легче с wavein, я был бы рад использовать это - если вы можете сказать мне, как).Затем я использую библиотеку FFTW - я знаю, как выполнить 1D, 2D (реальное и сложное) FFT, но я не совсем уверен, как это сделать, так как мне пришлось бы делать 3D FFT, чтобы получить частоту, амплитуду (это определитградиент цвета) и время.Или это просто 2D БПФ, и я получаю амплитуду и частоту?

Answer 1

Я использую скользящий ДПФ, который в много раз * в 1002 * раз быстрее, чем БПФ, в случае, когда вам нужно выполнять преобразование Фурье каждый раз, когда выборка поступает во входной буфер.

Это основано на том факте, что после выполнения преобразования Фурье для последних N выборок и получения новой выборки вы можете "отменить" эффект самой старой выборки и применить эффект последней выборки в один проход через данные Фурье! Это означает, что производительность скользящего ДПФ составляет O (n) по сравнению с O (Log2 (n) умножить на n) для БПФ. Кроме того, нет ограничений на степень двойки для размера буфера для поддержания производительности.

В полной программе тестирования, приведенной ниже, сравнивается скользящий ДПФ с fftw. В моем производственном коде я оптимизировал приведенный ниже код до нечитабельности, чтобы сделать его в три раза быстрее.

    #include <complex>
#include <iostream>
#include <time.h>
#include <math_defines.h>
#include <float.h>

#define DO_FFTW // libfftw
#define DO_SDFT

#if defined(DO_FFTW)
#pragma comment( lib, "d:\\projects\\common\\fftw\\libfftw3-3.lib" )

namespace fftw {
#include <fftw/fftw3.h>
}

fftw::fftw_plan plan_fwd;
fftw::fftw_plan plan_inv;

#endif

typedef std::complex<double> complex;

// Buffer size, make it a power of two if you want to improve fftw
const int N = 750;


// input signal
complex in[N];

// frequencies of input signal after ft
// Size increased by one because the optimized sdft code writes data to freqs[N]
complex freqs[N+1];

// output signal after inverse ft of freqs
complex out1[N];
complex out2[N];

// forward coeffs -2 PI e^iw -- normalized (divided by N)
complex coeffs[N];
// inverse coeffs 2 PI e^iw
complex icoeffs[N];

// global index for input and output signals
int idx;


// these are just there to optimize (get rid of index lookups in sdft)
complex oldest_data, newest_data;

//initilaize e-to-the-i-thetas for theta = 0..2PI in increments of 1/N
void init_coeffs()
{
    for (int i = 0; i < N; ++i) {
        double a = -2.0 * PI * i  / double(N);
        coeffs[i] = complex(cos(a)/* / N */, sin(a) /* / N */);
    }
    for (int i = 0; i < N; ++i) {
        double a = 2.0 * PI * i  / double(N);
        icoeffs[i] = complex(cos(a),sin(a));
    }
}


// initialize all data buffers
void init()
{
    // clear data
    for (int i = 0; i < N; ++i)
        in[i] = 0;
    // seed rand()
    srand(857);
    init_coeffs();
    oldest_data = newest_data = 0.0;
    idx = 0;
}

// simulating adding data to circular buffer
void add_data()
{
    oldest_data = in[idx];
    newest_data = in[idx] = complex(rand() / double(N));

}


// sliding dft
void sdft()
{
    complex delta = newest_data - oldest_data;
    int ci = 0;
    for (int i = 0; i < N; ++i) {
        freqs[i] += delta * coeffs[ci];
        if ((ci += idx) >= N)
            ci -= N;
    }
}

// sliding inverse dft
void isdft()
{
    complex delta = newest_data - oldest_data;
    int ci = 0;
    for (int i = 0; i < N; ++i) {
        freqs[i] += delta * icoeffs[ci];
        if ((ci += idx) >= N)
            ci -= N;
    }
}

// "textbook" slow dft, nested loops, O(N*N)
void ft()
{
    for (int i = 0; i < N; ++i) {
        freqs[i] = 0.0;
        for (int j = 0; j < N; ++j) {
            double a = -2.0 * PI * i * j / double(N);
            freqs[i] += in[j] * complex(cos(a),sin(a));
        }
    }
}

double mag(complex& c)
{
    return sqrt(c.real() * c.real() + c.imag() * c.imag());
}

void powr_spectrum(double *powr)
{
    for (int i = 0; i < N/2; ++i) {
        powr[i] = mag(freqs[i]);
    }

}

int main(int argc, char *argv[])
{
    const int NSAMPS = N*10;
    clock_t start, finish;

#if defined(DO_SDFT)
// ------------------------------ SDFT ---------------------------------------------
    init();

    start = clock();
    for (int i = 0; i < NSAMPS; ++i) {

        add_data();

        sdft();
        // Mess about with freqs[] here
        //isdft();

        if (++idx == N) idx = 0; // bump global index

        if ((i % 1000) == 0)
            std::cerr << i << " iters..." << '\r';
    }
    finish = clock();

    std::cout << "SDFT: " << NSAMPS / ((finish-start) / (double)CLOCKS_PER_SEC) << " fts per second." << std::endl;

    double powr1[N/2];
    powr_spectrum(powr1);
#endif

#if defined(DO_FFTW)

// ------------------------------ FFTW ---------------------------------------------
    plan_fwd = fftw::fftw_plan_dft_1d(N, (fftw::fftw_complex *)in, (fftw::fftw_complex *)freqs, FFTW_FORWARD, FFTW_MEASURE);
    plan_inv = fftw::fftw_plan_dft_1d(N, (fftw::fftw_complex *)freqs, (fftw::fftw_complex *)out2, FFTW_BACKWARD, FFTW_MEASURE);

    init();

    start = clock();
    for (int i = 0; i < NSAMPS; ++i) {

        add_data();

        fftw::fftw_execute(plan_fwd);
        // mess about with freqs here
        //fftw::fftw_execute(plan_inv);

        if (++idx == N) idx = 0; // bump global index

        if ((i % 1000) == 0)
            std::cerr << i << " iters..." << '\r';
    }
    // normalize fftw's output 
    for (int j = 0; j < N; ++j)
        out2[j] /= N;

    finish = clock();

    std::cout << "FFTW: " << NSAMPS / ((finish-start) / (double)CLOCKS_PER_SEC) << " fts per second." << std::endl;
    fftw::fftw_destroy_plan(plan_fwd);
    fftw::fftw_destroy_plan(plan_inv);

    double powr2[N/2];
    powr_spectrum(powr2);
#endif
#if defined(DO_SDFT) && defined(DO_FFTW)
// ------------------------------      ---------------------------------------------
    const double MAX_PERMISSIBLE_DIFF = 1e-11; // DBL_EPSILON;
    double diff;
    // check my ft gives same power spectrum as FFTW
    for (int i = 0; i < N/2; ++i)
        if ( (diff = abs(powr1[i] - powr2[i])) > MAX_PERMISSIBLE_DIFF)
            printf("Values differ by more than %g at index %d.  Diff = %g\n", MAX_PERMISSIBLE_DIFF, i, diff);

#endif
    return 0;
}

Answer 2

Если вам нужны амплитуда, частота и время на одном графике, то преобразование называется разложением по времени и частоте.Самый популярный из них называется «Кратковременное преобразование Фурье».Это работает следующим образом:
1. Возьмите небольшую часть сигнала (скажем, 1 секунду)
2. Окно с небольшим окном (скажем, 5 мс)
3. Вычислите 1-мерное преобразование Фурьеоконный сигнал.
4. Переместите окно на небольшую величину (2,5 мс)
5. Повторяйте вышеуказанные шаги до конца сигнала.
6. Все эти данные вводятся в матрицу, которая затем используетсясоздать вид трехмерного представления сигнала, который показывает его разложение по частоте, амплитуде и времени.

Длина окна определяет разрешение, которое вы можете получить в частотной и временной областях.Проверьте здесь для получения более подробной информации о STFT и найдите учебные пособия "Robi Polikar" по вейвлет-преобразованиям для ознакомления непрофессионала со всем вышеперечисленным.

Редактировать 1:
Вы берете функцию управления окнами (существует бесчисленное множество функций - вот список . Наиболее интуитивно понятным является прямоугольное окно, ночаще всего используются оконные функции Хэмминга / Хеннинга. Вы можете выполнить следующие шаги, если у вас в руке бумажный карандаш и нарисовать его вдоль.

Предположим, что полученный вами сигнал длится 1 секунду иназывается x[n]. Функция управления окном имеет длину 5 мс и называется w[n]. Поместите окно в начало сигнала (чтобы конец окна совпадал с точкой 5 мс сигнала) и умножьте x[n] и w[n] примерно так:
y[n] = x[n] * w[n] - умножение сигналов по точкам.
Возьмите БПФ y[n].

Затем вы немного сдвинете окно (скажем, 2,5 мс). Так что теперь окно растягивается от 2,5 мс до 7,5 мс сигнала x[n]. Повторите шаги умножения и генерации БПФ. Другими словами, у вас есть перекрытие 2,5 мсек. Вы увидите, что чаИзменение длины окна и перекрытия дает вам различные разрешения по оси времени и по частоте.

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

PS: Если вы поняли STFT и другие частотно-временные разложения сигнала, то у вас не должно быть проблем с шагами 2 и 4. То, что вы не поняли вышеупомянутые шаги, заставляет меня чувствовать, что выследует также пересмотреть частотно-временные разложения.

Answer 3

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

Answer 4

БПФ в реальном времени означает, что он полностью отличается от того, что вы только что описали.Это означает, что для данных N и X [N] ваш алгоритм дает Fx [i] при увеличении значения i .Это означает, что исходящее значение не будет вычислено, пока не будет завершено вычисление текущего значения.Это полностью отличается от того, что вы только что описали.

Аппаратное обеспечение обычно использует БПФ с количеством точек от 1 до 16 тысяч.Исправлено N, не вычисление в реальном времени.Перемещение окна БПФ, как описано в предыдущих ответах.