Member 13319218 Ответов: 1

Как выполнить БПФ с помощью библиотеки gsl на языке C


Дорогие друзья,

Я использую gsl для выполнения БПФ сложной функции с помощью тестового примера, приведенного в руководстве библиотеки gsl. Я успешно разработал свой код и его запуск, но не дал правильных результатов. Можете ли вы помочь мне решить эту проблему? Моя функция заключается в следующем:
double w = 0.1;
 double om = 0.7;
 double Q1 = 1;
 double Q2 = 2.1;
  double D1 = 2.1;

а функция, где я хочу применить БПФ, выглядит следующим образом
D=pow((cos(om*t1/2)),2)+((2*I*sin(om*t1)/om)*(Q2+pow(om,2)*Q1))-(16*D1*((pow(sin(om*t1/2),2))));


Что я уже пробовал:

#include <stdio.h>
#include <math.h>
#include <gsl gsl_errno.h>
#include <gsl gsl_fft_complex.h>
#include<complex.h>
#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])

int
main (void)
{
  int i;
 double w = 0.1;
 double om = 0.7;
 double Q1 = 1;
 double Q2 = 2.1;
int r = 1;
double D1 = 2.1;
  const int n = 630;
  double data[2*n];
  double c[2*n];
  double complex D; 
  gsl_fft_complex_wavetable * wavetable;
  gsl_fft_complex_workspace * workspace;

  for (i = 0; i < n; i++)
    {
      REAL(data,i) = 0.0;
      IMAG(data,i) = 0.0;
    }

  data[0] = 1.0;

  for (i = 1; i <= 500; i++)
    {
 double t1 = (m-1)*4.14;
 D=pow((cos(om*t1/2)),2)+((2*I*sin(om*t1)/om)*(Q2+pow(om,2)*Q1))-(16*D1*((pow(sin(om*t1/2),2))));
      REAL(data,i) = REAL(data,n-i) = creal(D); //r*cos(w*i);
      IMAG(data,i) = IMAG(data, n-i) = cimag(D); //r*sin(w*i);
    }

  for (i = 0; i < n; i++)
    {
      //printf ("%d: %e %e\n", i, REAL(data,i), 
//                                IMAG(data,i));
    }
  printf ("\n");

  wavetable = gsl_fft_complex_wavetable_alloc (n);
  workspace = gsl_fft_complex_workspace_alloc (n);

  for (i = 0; i < (int) wavetable->nf; i++)
    {
       printf ("# factor %d: %zu\n", i, 
               wavetable->factor[i]);
    }

  gsl_fft_complex_forward (data, 1.5, n, 
                           wavetable, workspace);

  for (i = 0; i < n; i++)
    {
     // printf ("%d: %e %e\n", i, REAL(data,i), 
//                                IMAG(data,i));
 c[i]  = REAL(data,i) + IMAG(data,i);    
printf("%e\n",c[i]);
}

  gsl_fft_complex_wavetable_free (wavetable);
  gsl_fft_complex_workspace_free (workspace);

   FILE *fs = fopen("spe.txt","w");

        for (i = 1; i<630; i++)
        {
        //double t = dx*(i-1);

        fprintf(fs," %d  %f\n", i, data[i]);
        }
        fclose(fs);
}

1 Ответов

Рейтинг:
0

KarstenK

Сравнивая ваш код с Научная библиотека GNU выглядит правильно, но они использовали фактор 1, а Вы 1,5 для int. Переосмыслите это.

Таким образом, ваш D выглядит немного странно, потому что это массив, но вы не обращаетесь к нему. Также результат правой части формулы выглядит постоянным. Может быть некоторые опечатка в переводе???

Итак, вот мое "предложение" для цикла:

for (i = 1; i <= 500; i++)
{
  // access the array
  D[i] = pow((cos(om*t1/2)),2)+((2*i/*use the loop counter*/*sin(om*t1)/om)*(Q2+pow(om,2)*Q1))-(16*D1*((pow(sin(om*t1/2),2))));
}


Member 13319218

Дорогой Карстенк,

Спасибо за Ваш быстрый ответ, но мой D уже является массивом t, так что нет необходимости делать его массивом "i".
И фактор не так уж сильно влияет,я его перепроверил.

Ваш самый радушный прием, если у вас есть еще какие-то предложения по этому поводу