Я хочу сделать БПФ аудиосигнала в реальном времени, то есть в то время как человек говорит в микрофоне. Я буду получать данные (я делаю это с помощью portaudio, если было бы проще с wavein, я был бы счастлив использовать это - если вы можете мне рассказать). Затем я использую библиотеку FFTW - я знаю, как выполнять 1D, 2D (real & complex) FFT, но я не уверен, как это сделать, так как мне нужно будет делать 3D FFT для получения частоты, амплитуды (это будет определять градиент цвета) и время. Или это просто 2D БПФ, и я получаю амплитуду и частоту?
Выполнение БПФ в режиме реального времени
Ответ 1
Если вам нужна амплитуда, частота и время на одном графике, тогда преобразование известно как разложение по частоте. Наиболее популярным называется короткосрочное преобразование Фурье. Он работает следующим образом:
1. Возьмите небольшую часть сигнала (скажем, 1 секунду)
2. Окончите его небольшим окном (скажем, 5 мс)
3. Вычислите 1D преобразование Фурье оконного сигнала.
4. Переместите окно на небольшое количество (2,5 мс)
5. Повторяйте выше шагов до конца сигнала.
6. Все эти данные вводятся в матрицу, которая затем используется для создания трехмерного представления сигнала, который показывает его разложение по частоте, амплитуде и времени.
Длина окна будет определять разрешение, которое вы можете получить в частотных и временных областях. Проверьте здесь для получения более подробной информации о STFT и найдите учебники "Robi Polikar" по вейвлет-преобразованиям для введения непрофессионала в вышесказанное.
Изменить 1:
Вы выполняете функцию окна (там есть неисчислимые функции - вот список. Наиболее интуитивно понятным является прямоугольное окно, но наиболее часто используемые окна Хэмминга/Ханнинга. Вы можете выполнить следующие шаги, если у вас есть бумажный карандаш в руке и нарисуйте его.
Предположим, что полученный вами сигнал длится 1 секунду и называется x[n]. Функция windowing длится 5 мс и называется w[n]. Поместите окно в начале сигнала (так что конец окна совпадает с точкой 5 мс сигнала) и умножьте x[n] и w[n] следующим образом: y[n] = x[n] * w[n] - умножение точек по точкам.
Возьмите БПФ y[n].
Затем вы меняете окно на небольшое количество (скажем, 2,5 мс). Итак, теперь окно простирается от 2,5 мс до 7,5 мс сигнала x[n]. Повторите шаги генерации умножения и FFT. Другими словами, у вас есть перекрытие 2,5 мс. Вы увидите, что изменение длины окна и перекрытия дает разные разрешения по времени и оси частот.
Как только вы это сделаете, вам нужно передать все данные в матрицу, а затем отобразить их. Наложение заключается в минимизации ошибок, которые могут возникнуть на границах, а также для получения более согласованных измерений за такие короткие временные рамки.
PS: Если вы понимали STFT и другие временные разложения сигнала, тогда у вас не должно было быть проблем с шагами 2 и 4. Чтобы вы не поняли вышеописанные шаги, мне кажется, что вы должны пересмотреть время -частотные разложения.
Ответ 2
Я использую раздвижную DFT, которая много раз быстрее, чем FFT в случае, когда вам нужно сделать преобразование Фурье каждый раз, когда образец поступает во входной буфер.
Это основано на том факте, что как только вы выполнили преобразование Фурье для последних N выборок, и придет новый образец, вы можете "отменить" эффект самого старого образца и применить эффект последнего образца в один проход через данные Фурье! Это означает, что скорость скольжения DFT равна O (n) по сравнению с O (Log2 (n) раз n) для FFT. Кроме того, нет ограничений на полномочия двух для размера буфера для поддержания производительности.
В полной программе испытаний ниже сравнивается скользящий DFT с 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 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;
}
Ответ 3
Вы можете создать БПФ в реальном времени, выбирая короткий промежуток времени и анализируя (FFT'ing) только этот период времени. Вероятно, вам удастся просто выбрать неперекрывающиеся временные интервалы, скажем, 100-500 миллисекунд; аналитически более чистый способ сделать это будет использовать скользящее окно (опять же, например, 100-500 мс), но это часто не нужно, и вы можете показывать красивую графику с неперекрывающимися временными интервалами без большой вычислительной мощности.
Ответ 4
В реальном времени FFT означает совершенно иное, чем то, что вы только что описали. Это означает, что при заданных N и X [N] ваш алгоритм дает Fx [i] при увеличении значения i. Значение, исходящее значение не вычисляется до завершения текущего вычисления значения. Это полностью отличается от того, что вы только что описали.
Аппаратное обеспечение обычно использует БПФ с отметками около 1 кБ. Исправлено N, а не вычисление в реальном времени. Перемещение окна FFT, как описано в предыдущих ответах.