如何在C中使用gsl库执行fft [英] How to perform fft using gsl library in C

问题描述

亲爱的朋友们，



我借助gsl库手册中给出的测试示例，使用gsl执行复杂函数的fft。我已经成功开发了我的代码并且它正在运行但没有产生正确的结果。你能帮帮我解决这个问题吗？我的功能如下：

Dear folks,

I am using gsl to perform fft of a complex function with help of a test example given in manual of gsl library. I have developed my code successfully and its running but not producing correct results. Can you help me out to resolve this problem? My function is as follows:

double w = 0.1;
 double om = 0.7;
 double Q1 = 1;
 double Q2 = 2.1;
  double D1 = 2.1;

和函数，我想申请fft，如下：

and function, where I want to apply fft, is as follow

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))));

我的尝试：

What I have tried:

#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);
}

推荐答案

将您的代码与GNU科学图书馆看起来是正确的，但他们使用了因子1，但你使用的是1.5。重新考虑一下。



所以你的D看起来有点奇怪，因为它是一个数组，但是没有访问它。此外，公式右侧的结果看起来是不变的。也许在翻译中有一些拼写错误 ???



所以这是循环的建议：

By comparing your code with the GNU scientific library is looking correct, but they have used the factor 1 but you 1.5 for an int. Rethink that.

So your D is looking a bit strange, because it is an array, but arent accessing it. Also the result of the right side of the formula looks constant. Maybe some typo in the translation???

So here is my "proposal" for the loop:

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))));
}

这篇关于如何在C中使用gsl库执行fft的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！