Пытаюсь спектрограмму сигнала нарисовать, есть вот такой тестовый код:
#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:-Вот эта функция
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])));
}
}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);
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);
}
