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

#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>
#include <math.h>
#include <complex.h>
#include <fftw3.h>

#define SQR(x) ((x)*(x))
#define XRESOL 1024
#define YRESOL 768

uint8_t pic[YRESOL][XRESOL] = {0};

#define FFT_CHUNK 512

void fft_calculate(double in[FFT_CHUNK], double complex out[(FFT_CHUNK/2)+1])
{
  fftw_plan p = fftw_plan_dft_r2c_1d(FFT_CHUNK, in, out, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
  fftw_execute(p);
  fftw_destroy_plan(p);
}

void convert(double complex in[(FFT_CHUNK/2)+1], double out[(FFT_CHUNK/2)+1])
{
  for (size_t i = 0; i < (FFT_CHUNK/2)+1; i++)
  {
    out[i] = sqrt(SQR(creal(in[i])) + SQR(cimag(in[i])));
  }
}

int main(void)
{
  double arr[XRESOL];
  for (size_t i = 0; i < XRESOL; i++)
  {
    arr[i] =
      sin((double)i * M_2_PI * 0.2)*0.1 +
      sin((double)i * M_2_PI * 0.4)*0.1 +
      sin((double)i * M_2_PI * 0.6)*0.1 +
      sin((double)i * M_2_PI * 0.8)*0.1 +
      sin((double)i * M_2_PI * 1.0)*0.1;
  }
  for (size_t i = 0; i < XRESOL - FFT_CHUNK; i++)
  {
    double complex tmp[(FFT_CHUNK/2)+1];
    double result[(FFT_CHUNK/2)+1];
    fft_calculate(&arr[i], tmp);
    convert(tmp, result);
    for (size_t j = 0; j < (FFT_CHUNK/2)+1; j++)
    {
      #define MULT 10
      pic[j][i] = (result[j] * MULT) > 0xff ? 0xff : (uint8_t)(result[j]* MULT);
      #undef MULT
    }
  }

  fwrite(pic, sizeof(pic), 1, stdout);
  return EXIT_SUCCESS;
}

Собрать-запустить можно вот так:
gcc main.c -lfftw3 -lm
./a.out | display -size 1024x768 -depth 8 gray:-

И оно нарисует три полоски.
Понятно, что в функции fft_calculate() можно переиспользовать план, как указано в этих доках http://www.fftw.org/fftw3_doc/New_002darray-Execute-Functions.html а не создавать каждый раз. Но вот как там с real to real функциями из этого fftw работать? И нужно ли их вообще тут использовать?

Вот эта функция

void convert(double complex in[(FFT_CHUNK/2)+1], double out[(FFT_CHUNK/2)+1])
{
  for (size_t i = 0; i < (FFT_CHUNK/2)+1; i++)
  {
    out[i] = sqrt(SQR(creal(in[i])) + SQR(cimag(in[i])));
  }
}

Правильно ли это? Насколько это вообще корректно? Может тут стоит использовать real to real функции, а не суммировать квадраты действительной и мнимой части, извлпкая потом квадратный корень?

http://www.fftw.org/fftw3_doc/Real_002dto_002dReal-Transforms.html для real to real преобразования там есть параметр fftw_r2r_kind

fftw_plan fftw_plan_r2r_1d(int n, double *in, double *out,
                           fftw_r2r_kind kind, unsigned flags);

Там есть много этих кайндов: http://www.fftw.org/fftw3_doc/Real_002dto_002dReal-Transform-Kinds.html какой из них мне надо использовать? Я попробовал так сделать
void fft_calculate(double in[FFT_CHUNK], double out[FFT_CHUNK])
{
  fftw_plan p = fftw_plan_r2r_1d(FFT_CHUNK, in, out, FFTW_DHT, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);
  fftw_execute(p);
  fftw_destroy_plan(p);
}

но получилась какая-то ерунда








 , , , ,






URL записи