用小波变换对数据多尺度分析，写了个程序可以正常运行，但运行结果好像有错误（每个尺度下的小波系数的前几个数和后几个数特别大），不知是何原因，望大神指教。
#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 编辑 ]
