Векторы и матрицы в 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 блоков, в результате чего получится двумерный массив:

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), обычно используемая для преобразования линейной шкалы в децибелы, использует десятичный логарифм вместо естественная функция 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;

Похоже, вы вообще пропустили начальное объявление для «v», а «in» не объявлено должным образом.

См. эту страницу для связанного вопрос о создании 2D-массивов в C++. Насколько я понимаю, fftw_malloc() в основном является new() или malloc(), но правильно выравнивает переменную для алгоритма FFTW.

Поскольку вы не предоставляете 'v' чему-либо, связанному с FFTW, вы можете использовать для этого стандартный malloc().