傅里叶级数系数​​ gnuplot

Posted 2023-03-31

【中文标题】傅里叶级数系数​​ gnuplot

【英文标题】：Fourier series coefficients gnuplot

【发布时间】：2015-12-02 12:19:41

【问题描述】：

这个程序应该找到傅立叶级数的系数，然后绘制从计算值构建的函数。

分段是

h(t)=1+(t/pi) 对于 t0

#include <stdio.h>
#include <math.h>
#include "comphys.c"
#include "comphys.h"
#define pi 3.141592653589793

double trap1 (double l, double u, int m, int N);
double trap(double l,double u,int m,int N); // Use trapezoidal method
double inte(double x,int m);
double integrand(double x, int m); // This is the integrand of the bms
double w=1.0; //omega (here - w-2pi/Period where Period=2pi)
int main()

  FILE *out,*rsp;
  int m; // Index of the b coefficients
  int M=5; // Number of coeeficients to compute
  int T=201;  // Number of time steps 
  int t; // This is a counter to keep track of the h(t) function to plot
  double *b,*h,j,dj,*a,a0;
  double ll,ul;
  a=dvector(1,M);
  b=dvector(1,M);  // coefficients
  h=dvector(1,T);  // function built from Fourier terms
  ll=-pi;  // lower time limit
  ul=pi;   // upper time limit

  // Prepare to write for plotting
  printf("\nBrute force complex Fourier transform algorithm");
  if ((out=fopen("fs.out","w"))==NULL) 
    printf("\nCannot open file for output\n");
    getchar();
    return(0);          
  

  // Loop through the wave numbers    
  for (m=1;m<=M;m++)
    
      b[m]=(1.0/pi)*trap(ll,ul,m,T);
      a[m]=(1.0/pi)*trap1(ll,ul,m,T);
      printf("\nb[%d]=%f   a[%d]=%f",m,b[m],m,a[m]);

      // Build the h(t) function as we go
      t=0;
      dj=(ul-ll)/(double)T;
      for (j=ll;j<=ul;j+=dj)
    
      t++;
      h[t]+=(b[m]*sin((double)m*w*j))+(a[m]*cos((double)m*w*j));
    
    

  // Write out for plotting
  t=0;
  for (j=ll;j<=ul;j+=dj)
    
      t++;
      fprintf(out,"%f %f\n",j,h[t]);
    
  fclose(out);
  printf("\nProgram complete without known error.\n");

  if((rsp = fopen("gnuplot.rsp","w")) == NULL)
    
      printf("\nCannot open gnuplot.rsp for writing\n");
      exit(1);
    
  fprintf(rsp,"plot '%s' using 1:2 with lines\n",out);
  fprintf(rsp,"pause mouse\n");
  fprintf(rsp,"replot\n");
  fclose(rsp);
  if(system("gnuplot gnuplot.rsp") == -1)
    
      printf("\nCommand could not be executed\n");
      exit(1);
    
  return;


double trap (double l, double u, int m, int N)


  double x,integral,dx;
  int i;

  if (u==l) return (0.0);
  dx=(u-l)/(double)(N-1);

  integral=0.0;
  integral+=0.5*(integrand(u,m)-integrand(l,m));
  x=l+dx;
  for (i=1;i<N;i++)
    
      integral += integrand(x,m);
      x+=dx;
    
  integral*=dx;
  printf("\n%f\n",integral);
  return(integral);


double trap1 (double l, double u, int m, int N)

  double x,integral,dx;
  int i;

  if (u==l) return (0.0);
  dx=(u-l)/(double)(N-1);

  integral=0.0;
  integral+=0.5*(inte(u,m)-inte(l,m));
  x=l+dx;
  for (i=1;i<N;i++)
    
      integral += inte(x,m);
      x+=dx;
    
  integral*=dx;
  return(integral);



double integrand(double x, int m)

  if (x<0) return (sin((double)m*w*x)+((x/pi)*sin((double)m*w*x)));
  return(sin((double)m*w*x)-((x/pi)*sin((double)m*w*x)));


double inte(double x,int m)

  if (x<0) return (cos((double)m*w*x)+((x/pi)*cos((double)m*w*x)));
  return(cos((double)m*w*x)-((x/pi)*cos((double)m*w*x)));

我遇到的问题是，它编译得很好并列出了系数（不知道它们是否正确或任何东西），但最后在尝试绘制时它给了我

"无法读取数据文件"?,??" "gnuplot.rsp", line:1 util.c: 没有这样的文件或目录"

我不知道如何解决这个问题。我用 gnuplot 绘制了很多程序，所以它不应该是丢失的文件或任何东西。

【问题讨论】：

你确定你在gnuplot.rsp所在的同一目录中运行程序吗？看来您根本没有包含名为 so 的文件的当前工作目录。

【参考方案1】：

这里：

fprintf(rsp,"plot '%s' using 1:2 with lines\n",out);  

你应该直接使用字符串“fs.out”而不是“out”，它是一个文件指针（也已经关闭）。