fft — Векторы и матрицы в C ++ для генерации спектрограммы

Как указано в Документация FFTW, размер вывода (out в вашем случае) при использовании fftw_plan_dft_r2c_1d не совпадает с размером ввода. Более конкретно для ввода N реальные образцы, выход состоит из N/2+1 сложные ценности. Затем вы можете выделить out с:

out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2 + 1));

Для вывода спектрограммы у вас будет аналогично (N/2+1) величины для каждого из B блоки, в результате чего в 2D массив:

double** v = new double*[B];
for (int i = 0; i < B; i++){
v[i] = new double[(N/2+1)];
}

Также обратите внимание, что вы можете повторно использовать входной буфер in для каждой итерации (заполняя ее данными для нового блока). Однако, так как вы решили вычислить NТочка БПФ и будет хранить меньшие блоки Windowsize образцы (в этом случае N=500 а также Windowsize=100), не забудьте инициализировать оставшиеся образцы нулями:

in = (double*)fftw_malloc(sizeof(double) * N);
for (int i = 0; i < N; i++){
in[i] = 0;
}

Обратите внимание, что в дополнение к объявлению и распределению in а также v переменных, код, который вы разместили, страдает от нескольких дополнительных проблем:

• При вычислении Окно Ханнинга, вы должны разделить на Windowsize-1 не N-1 (так как в вашем случае N соответствуют размеру БПФ).
• Вы берете БПФ того же блока signal снова и снова, так как вы всегда индексируете j в [0,Windowsize] спектр. Скорее всего, вы захотите добавить смещение каждый раз, когда обрабатываете другой блок.
• Поскольку размер БПФ не изменяется, вам нужно создать план только один раз. По крайней мере, если вы собираетесь создавать свой план на каждой итерации, вы должны аналогичным образом уничтожить его (с fftw_destroy_plan) на каждой итерации.

И несколько дополнительных моментов, которые могут потребовать некоторых мыслей:

• Масштабирование масштабируемых по логарифму величин путем деления на N может не делать то, что вы думаете. Вы, скорее всего, захотите масштабировать линейные масштабы (т.е. разделите величину, прежде чем брать логарифм). Обратите внимание, что это приведет к постоянному смещению кривой спектра, что для многих приложений не столь существенно. Если масштабирование важно для вашего приложения, вы можете взглянуть на другой мой ответ Больше подробностей.
• Общая формула 20*log10(x) обычно используется для преобразования линейного масштаба в децибел использует логарифм основания-10 вместо естественного log (база e~2.7182) функция, которую вы использовали. Это приведет к мультипликативному масштабированию (растяжению), которое может или не может быть значительным в зависимости от вашего приложения.

Подводя итог, следующий код может более соответствовать тому, что вы пытаетесь сделать:

// Allocate & initialize buffers
in = (double*)fftw_malloc(sizeof(double) * N);
for (int i = 0; i < N; i++){
in[i] = 0;
}
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2 + 1));
v = new (double*)[B];
for (int i = 0; i < B; i++){
v[i] = new double[(N/2+1)];
}

// Generate the signal
...

// Create the plan once
plan_forward = fftw_plan_dft_r2c_1d (Windowsize, in, out, FFTW_ESTIMATE);

// Applying the Hanning window function on each block B
for(int k = 0; k < B; k++){
for(int j = 0; j < Windowsize; j++){
// Hanning Window
double multiplier = 0.5 * (1 - cos(2 * M_PI * j / (Windowsize-1)));
in[j] = multiplier * signal[j+k*Windowsize];
}
fftw_execute(plan_forward);
for (int j = 0; j <= N/2; j++){
// Factor of 2 is to account for the fact that we are only getting half
// the spectrum (the other half is not return by a R2C plan due to symmetry)
v[k][j] = 2*(out[j][0] * out[j][0] + out[j][1] * out[j][1])/(N*N);
}
// DC component and at Nyquist frequency do not have a corresponding symmetric
// value, so should not have been doubled up above. Correct those special cases.
v[k][0]   *= 0.5;
v[k][N/2] *= 0.5;
// Convert to decibels
for (int j = 0; j <= N/2; j++){
// 20*log10(sqrt(x)) is equivalent to 10*log10(x)
// also use some small epsilon (e.g. 1e-5) to avoid taking the log of 0
v[k][j] = 10 * log10(v[k][j] + epsilon);
}
}

// Clean up
fftw_destroy_plan(plan_forward);
fftw_free(in);
fftw_free(out);
// Delete this last one after you've done something useful with the spectrogram
for (int i = 0; i < B; i++){
delete[] v[i];
}
delete[] v;