#include "stdio.h"    //标准输入（scanf）输入（printf）头文件
  #include "string.h"   //字符数组的函数定义的头文件，常用函数有strlen、strcmp、strcpy等等
  #include "malloc.h"   //动态存储分配函数头文件,当对内存区进行操作时,调用相关函数
  #include "stdlib.h"   //标准函数库的头文件，声明了数值与字符串转换函数，伪随机数生成函数，动态内存分配函数，进程控制函数等公共函数
  #include "math.h"     //数学函数库的头文件,一些数学计算的公式的具体实现是放在math.h里
  #include "conio.h"    //定义了通过控制台进行数据输入和数据输出的函数，主要是一些用户通过按键盘产生的对应操作，比如getch()函数等等

  #define channel_number 19        //导联数
  #define sample_number 20480      //样本数据
  #define person_number 6          //文件数
  #define scale_number  30         //尺度数
  #define first_step 1
  #define first_scale 0.002        //决定了可辨识的最大频率为100
  #define scale_coefficient 1.2    //迭代步长，和（scale_number、first_scale）一起决定了可辨识的最小频率为0.5
  #define sample_time 20.0         //采样时间

  void stradd(char *ch,int jj);

  char fn1[20]={"np"},    //np是data文件名
       fn2[20]={"ws"},    //ws是子波系数文件名
       fn3[20]={"wp"},   
       fn4[20]={"signal"},
       fn5[20]={"report1"},  //能量文件名
       fn6[20]={"report2"};

  FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp6;

  double u[20480][20],t1=0.0009765625,a[31],es[20][31],*c,*g,th,b,s,e,sum_new,sum_old;

  main()
  {
      int i,j,k,k1,k2,k3,k4,ejection_num,sweep_number,search_step;

      if((c=(double *)malloc(20480*sizeof(double)))==NULL)
      {
          printf( " malloc can not be opened ");
          getch();
          exit(0);
      }

      if((g=(double *)malloc(20480*sizeof(double)))==NULL)
      {
          printf( " malloc can not be opened ");
          getch();
          exit(0);
      }

      printf("欢迎您进入Gauss子波分析程序\n");

      //t1=0.0009765625;    //采样周期，即采样频率的倒数
      a[1]=first_scale;
      search_step=first_step;

      stradd(fn1,0);
      stradd(fn2,0);
      stradd(fn3,0);
      stradd(fn4,0);
      stradd(fn5,0);
      stradd(fn6,0);

      for(i=1;i<=person_number;i++)
      {
          printf("The person number is %d\n",i);

          stradd(fn1,i);
          stradd(fn2,i);
          stradd(fn3,i);
          stradd(fn4,i);
          stradd(fn5,i);
          stradd(fn6,i);

          if((fp1=fopen(fn1,"r"))==NULL)
          {
              printf("Can not open fn1!\n");
              getch();
              exit(0);
          }
   
          if((fp4=fopen(fn4,"w"))==NULL)
          {
              printf("Can not open fn4!\n");
              getch();
              exit(0);
          }
   
          if((fp5=fopen(fn5,"w"))==NULL)
          {
              printf("can not open fn5!\n");
              getch();
              exit(0);
          }
   
          if((fp6=fopen(fn6,"w"))==NULL)
          {
              printf("can not open fn6!\n");
              getch();
              exit(0);
          }

          for(j=0;j<sample_number;j++)
          {
              fprintf(fp4,"%.5f\t",j*t1);
              for(k=1;k<=channel_number;k++)
              {
                  fscanf(fp1,"%lf",&u[j][k]);
                  fprintf(fp4,"%.1f\t",u[j][k]);
              }
              fprintf(fp4,"\n");
          }
          fclose(fp1);
          fclose(fp4);

          stradd(fn2,0);
          stradd(fn3,0);
         
          e=0.0;

          for(k1=1;k1<=channel_number;k1++)
          {
              printf("  The channel number is %d\n",k1);

              stradd(fn2,k1);
              stradd(fn3,k1);

              if((fp3=fopen(fn3,"w"))==NULL)
              {
                  printf("can not open fn3!\n");
                  getch();
                  exit(0);
              }

              for(k3=0;k3<sample_number;k3++)
              {
                  c[k3]=u[k3][k1];
              }
              stradd(fn2,0);
   
              for (k2=1;k2<=scale_number;k2++)
              {
                  stradd(fn2,k2);
     
                  if((fp2=fopen(fn2,"w"))==NULL)
                  {
                      printf("can not open fn2!\n");
                      getch();
                      exit(0);
                  }

                  th=t1/a[k2];
                  es[k1][k2]=0.0;
                  b=0.0;
                  s=0.0;
                  ejection_num=0;
                  sweep_number=0;
                  sum_old=0.0;
                  
                  for (k3=0;k3<sample_number;k3++)
                  {
                      s=k3*th;
                      g[k3]=sqrt(2)*s*exp(-0.5*s*s)/sqrt(sqrt(3.14));
                      //s+=th;
                      //printf("%d\t%e\n",k3,g[k3]);
                  }
                  
                  for (k3=0;k3<sample_number;k3+=search_step)
                  {
                      sum_new=0.0;

                      for (k4=0;k4<k3;k4++)
                      {
                          sum_new+=-c[k4]*g[k3-k4];
                          sum_new+=c[k4]*g[sample_number-k3+k4];
                      }
      
                      for (k4=k3;k4<sample_number;k4++)
                      {
                          sum_new+=c[k4]*g[k4-k3];
                          sum_new+=-c[k4]*g[sample_number+k3-k4-1];
                      }
                      sum_new=sum_new/sqrt(a[k2]);
                      //printf("%e\n",sum_new);
      
                      if((sum_old>0)&&(sum_new<0))
                      {
                          ejection_num+=1;
                      }
      
                      if((sum_old<0)&&(sum_new>0))
                      {
                          sweep_number+=1;
                      }

                      fprintf(fp3,"%d\t%d\t%f\n",k3,k2,sum_new);
                      fprintf(fp2,"%d\t%f\n",k3,sum_new);
                      //printf("%d\t%e\n",k3,sum_new);


                     es[k1][k2]+=sum_new*sum_new;
                     b=b+t1*search_step;
                     sum_old=sum_new;
                  }
                  fclose(fp2);
                  a[k2+1]=a[k2]*scale_coefficient;
                  e+=es[k1][k2];
                  printf("    The scale number is %d\t%f\t%f\n",k2,a[k2],es[k1][k2]);
              }
              fclose(fp3);

              stradd(fn2,-1);
          }
          stradd(fn2,-1);
          stradd(fn3,-1);

          for (k2=1;k2<=scale_number;k2++)
          {
              fprintf(fp5,"%d\t%e\t%e\t",k2,a[k2],0.2/a[k2]);

              for(k1=1;k1<=channel_number;k1++)
              {
                  es[k1][k2]=es[k1][k2]/e;
                  fprintf(fp5,"%e\t",es[k1][k2]);
                  fprintf(fp6,"%d\t%d\t%e\n",k1,k2,es[k1][k2]);
              }
              fprintf(fp5,"\n");
          }
          fclose(fp5);
          fclose(fp6);
      }
      free(c);
      free(g);

      printf("The program is coming to the end!");
      getch();
      return 0;
  }


  void stradd(char *ch,int i)
  {     
      int i1,i2,n;

      n=strlen(ch);

      switch(i)
      {     
          case -1 :   
              ch[n-1]=0;
                 ch[n-2]=0;   
          break;
      
          case 0 :   
              ch[n]='0';
                 ch[n+1]='0';   
          break;

          default:   
              i1=i/10+'0';         // i/10：表示i整除10，i%10:表示i除以10后取余数。
                i2=i-i/10*10+'0';
              ch[n-2]=(char)i1;
              ch[n-1]=(char)i2;      
          break;
      }
  }

PS：有6个数据文件，每个文件有19列，每列为一个信号，采样频率是1024Hz，时长20s。

[ 本帖最后由 driverfeng 于 2015-7-29 11:43 编辑 ]