小波去噪c语言程序

嵌入式设计应用

127人已加入

描述

1、小波阈值去噪理论

小波阈值去噪就是对信号进行分解,然后对分解后的系数进行阈值处理,最后重构得到去噪信号。该算法其主要理论依据是:小波变换具有很强的去数据相关性,它能够使信号的能量在小波域集中在一些大的小波系数中;而噪声的能量却分布于整个小波域内。因此,经小波分解后,信号的小波系数幅值要大于噪声的系数幅值。可以认为,幅值比较大的小波系数一般以信号为主,而幅值比较小的系数在很大程度上是噪声。于是,采用阈值的办法可以把信号系数保留,而使大部分噪声系数减小至零。小波阈值收缩法去噪的具体处理过程为:将含噪信号在各尺度上进行小波分解,设定一个阈值,幅值低于该阈值的小波系数置为0,高于该阈值的小波系数或者完全保留,或者做相应的“收缩(shrinkage)”处理。最后将处理后获得的小波系数用逆小波变换进行重构,得到去噪后的信号.

2、小波阈值去噪c语言程序

此程序是用于信号处理分析,突出奇异值的前段处理,对信号进行小波包分解,用C语言实现的,仅供参考。

#include

#include

#include

#include

#defineLENGTH4096//信号长度

#defineDB_LENGTH8//Daubechies小波基紧支集长度

/******************************************************************

* 一维卷积函数

* 说明:循环卷积,卷积结果的长度与输入信号的长度相同

* 输入参数:data[],输入信号;core[],卷积核;cov[],卷积结果;

*n,输入信号长度;m,卷积核长度。

******************************************************************/

/*voidCovlution(doubledata[],doublecore[],doublecov[],intn,intm)

{

inti=0;

intj=0;

intk=0;

//将cov[]清零

for(i=0;i《n;i++)

{

cov[i]=0;

}

//前m/2+1行

i=0;

for(j=0;j《m/2;j++,i++)

{

for(k=m/2-j;k《m;k++)

{

cov[i]+=data[k-(m/2-j)]*core[k];//k针对core[k]

}

for(k=n-m/2+j;k《n;k++)

{

cov[i]+=data[k]*core[k-(n-m/2+j)];//k针对data[k]

}

}

//中间的n-m行

for(i=m/2;i《=(n-m)+m/2;i++)

{

for(j=0;j《m;j++)

{

cov[i]+=data[i-m/2+j]*core[j];

}

}

//最后m/2-1行

i=(n-m)+m/2+1;

for(j=1;j《m/2;j++,i++)

{

for(k=0;k《j;k++)

{

cov[i]+=data[k]*core[m-j-k];//k针对data[k]

}

for(k=0;k《m-j;k++)

{

cov[i]+=core[k]*data[n-(m-j)+k];//k针对core[k]

}

}

}

*/

//定义一个线性卷积

voidCovlution(doubledata[],doublecore[],doublecov[],intn,intm)

{

inti=0;

intj=0;

intt=0;

//将cov[]清零

for(j=0;j《n+m-1;j++)

{

cov[j]=0;

}

for(j=0;j《m+n-1;j++)

{

if(j《=m-1)//前面m行

{

for(i=0,t=j;t》=0;i++,t--)

cov[j]+=data[i]*core[t];

}

elseif(j《=n-1)//中间n-m行

{

for(i=j-m+1,t=m-1;t》=0;i++,t--)

cov[j]+=data[i]*core[t];

}

else//后面m行

{

for(i=j-m+1,t=m-1;i《n;i++,t--)

cov[j]+=data[i]*core[t];

}

}

}

/******************************************************************

* 一维小波变换函数

* 说明:一维小波变换,只变换一次

* 输入参数:input[],输入信号;output[],小波变换结果,包括尺度系数和

* 小波系数两部分;temp[],存放中间结果;h[],Daubechies小波基低通滤波器系数;

* g[],Daubechies小波基高通滤波器系数;n,输入信号长度;m,Daubechies小波基紧支集长度。

******************************************************************/

voidDWT1D_1(doubleinput[],doubleoutput0[],doubleoutput1[],doubletemp[],doubleh[],

doubleg[],intn,intm)

{

//doubletemp[LENGTH]={0};//?????????????

inti=0;

/*

//尺度系数和小波系数放在一起

Covlution(input,h,temp,n,m);

for(i=0;i《n;i+=2)

{

output[i]=temp[i];

}

Covlution(input,g,temp,n,m);

for(i=1;i《n;i+=2)

{

output[i]=temp[i];

}

*/

//尺度系数和小波系数分开

Covlution(input,h,temp,n,m);

for(i=0;i《n+m-1;i+=2)

{

output0[i/2]=temp[i];//尺度系数,进行了2抽值,即尺度空间,低频概貌部分

}

Covlution(input,g,temp,n,m);

for(i=1;i《n+m-1;i+=2)

{

//output[i]=temp[i];

output1[(i-1)/2]=temp[i];//小波系数,已经进行了2抽取,即高频细节部分

}

}

voidDWT1D_2(doubleinput[],doubleAA2[],doubleDA2[],doubleAD2[],doubleDD2[],doubletemp0[],doubletemp1[],doubletemp[],doubleh[],

doubleg[],intn,intm)

{

DWT1D_1(input,temp0,temp1,temp,h,g,n,m);

inta1=(m+n)/2;

DWT1D_1(temp0,AA2,DA2,temp,h,g,a1,m);

intd1=(n+m-4)/2;

DWT1D_1(temp1,AD2,DD2,temp,h,g,d1,m);

}

voidDWT1D_3(doubleinput[],doubleAAA3[],doubleDAA3[],doubleADA3[],doubleDDA3[],doubleAAD3[],doubleDAD3[],doubleADD3[],doubleDDD3[],

doubletemp0[],doubletemp1[],doubletemp2[],doubletemp3[],doubletemp00[],doubletemp11[],doubletemp[],doubleh[],doubleg[],intn,intm)

{

DWT1D_2(input,temp0,temp1,temp2,temp3,temp00,temp11,temp,h,g,n,m);

intaa2=(m+n)/4;

DWT1D_1(temp0,AAA3,DAA3,temp,h,g,aa2,m);

intda2=(n+3*m-8)/4;

DWT1D_1(temp1,ADA3,DDA3,temp,h,g,da2,m);

intad2=(n+3*m-4)/4;

DWT1D_1(temp2,AAD3,DAD3,temp,h,g,ad2,m);

intdd2=(n+3*m-12)/4;

DWT1D_1(temp3,ADD3,DDD3,temp,h,g,dd2,m);

}

voidmain()

{

doubledata[LENGTH];//输入信号

doubletemp0[(LENGTH+DB_LENGTH)/4];//先定义了一个中间结果

doubletemp1[(LENGTH+DB_LENGTH*3-8)/4];

doubletemp2[(LENGTH+DB_LENGTH*3-4)/4];

doubletemp3[(LENGTH+DB_LENGTH*3-12)/4];

doubletemp00[(LENGTH+DB_LENGTH)/2];

doubletemp11[(LENGTH+DB_LENGTH-4)/2];

doubletemp[LENGTH+DB_LENGTH-1];

doubleAAA3[(LENGTH+DB_LENGTH*5)/8];//一维小波变换后的结果

doubleDAA3[(LENGTH+DB_LENGTH*5-16)/8];

doubleADA3[(LENGTH+DB_LENGTH*7-8)/8];

doubleDDA3[(LENGTH+DB_LENGTH*7-24)/8];

doubleAAD3[(LENGTH+DB_LENGTH*7-4)/8];

doubleDAD3[(LENGTH+DB_LENGTH*7-20)/8];

doubleADD3[(LENGTH+DB_LENGTH*7-12)/8];

doubleDDD3[(LENGTH+DB_LENGTH*7-28)/8];

intaaa3=(LENGTH+DB_LENGTH*5)/8;//一维小波变换后的结果数组的长度

intdaa3=(LENGTH+DB_LENGTH*5-16)/8;

intada3=(LENGTH+DB_LENGTH*7-8)/8;

intdda3=(LENGTH+DB_LENGTH*7-24)/8;

intaad3=(LENGTH+DB_LENGTH*7-4)/8;

intdad3=(LENGTH+DB_LENGTH*7-20)/8;

intadd3=(LENGTH+DB_LENGTH*7-12)/8;

intddd3=(LENGTH+DB_LENGTH*7-28)/8;

intn=0;//输入信号长度

intm=8;//Daubechies正交小波基长度

inti=0;

chars[32];//从txt文件中读取一行数据

/*//DB3

staticdoubleh[]={.332670552950,.806891509311,.459877502118,-.135011020010,

-.085441273882,.035226291882};

staticdoubleg[]={.035226291882,.085441273882,-.135011020010,-.459877502118,

.806891509311,-.332670552950};

*/

//DB4

staticdoubleh[]={0.2303778133088964,0.7148465705529154,0.6308807679398587,-0.0279837694168599,-0.1870348117190931,0.0308413818355607,0.0328830116668852,-0.0105974017850690};//h[],Daubechies小波基低通滤波器系数;

staticdoubleg[]={-0.0105974017850690,-0.0328830116668852,0.0308413818355607,0.1870348117190931,-0.0279837694168599,-0.6308807679398587,0.7148465705529154,-0.2303778133088964};//g[],Daubechies小波基高通滤波器系数

//读取输入信号

FILE*fp;

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

fp=fopen(“heb1.txt”,“r”);

if(fp==NULL)//如果读取失败

{

printf(“错误!找不到要读取的文件hea1.txt/n”);

exit(1);//中止程序

}

while(fgets(s,32,fp)!=NULL)//读取长度n要设置得长一点,要保证读到回车符,这样指针才会定位到下一行?回车符返回的是零值?是,非数字字符经过atoi变换都应该返回零值

{

//fscanf(fp,“%d”,&data[count]);//一定要有“&”啊!!!最后读了个回车符!适应能力不如atoi啊

data[n]=atof(s);

n++;

}

//一维小波变换

//DWT1D(data,data_output,temp,h,g,n,m);

DWT1D_3(data,AAA3,DAA3,ADA3,DDA3,AAD3,DAD3,ADD3,DDD3,

temp0,temp1,temp2,temp3,temp00,temp11,temp,h,g,n,m);

//一维小波变换后的结果写入txt文件

fp0=fopen(“data_output_db4_aaa3.txt”,“w”);

fp1=fopen(“data_output_db4_daa3.txt”,“w”);

fp2=fopen(“data_output_db4_ada3.txt”,“w”);

fp3=fopen(“data_output_db4_dda3.txt”,“w”);

fp4=fopen(“data_output_db4_aad3.txt”,“w”);

fp5=fopen(“data_output_db4_dad3.txt”,“w”);

fp6=fopen(“data_output_db4_add3.txt”,“w”);

fp7=fopen(“data_output_db4_ddd3.txt”,“w”);

//打印一维小波变换后的结果

for(i=0;i《aaa3;i++)

{//if(i==0)

printf(“%s\n”,“AAA3”);

printf(“%f\n”,AAA3[i]);

fprintf(fp0,“%f\n”,AAA3[i]);

//}

//else

//{

//printf(“%f\n”,AAA3[i]);

//fprintf(fp0,“%f\n”,AAA3[i]);

//}

}

for(i=0;i《daa3;i++)

{//if(i==0)

printf(“%s\n”,“DAA3”);

printf(“%f\n”,DAA3[i]);

fprintf(fp1,“%f\n”,DAA3[i]);

//}

//else

//{

//printf(“%f\n”,DAA3[i]);

//fprintf(fp1,“%f\n”,DAA3[i]);

//}

}

for(i=0;i《ada3;i++)

{//if(i==0)

printf(“%s\n”,“ADA3”);

printf(“%f\n”,ADA3[i]);

fprintf(fp2,“%f\n”,ADA3[i]);

//}

//else

//{

//printf(“%f\n”,ADA3[i]);

//fprintf(fp2,“%f\n”,ADA3[i]);

//}

}

for(i=0;i《dda3;i++)

{//if(i==0)

printf(“%s\n”,“DDA3”);

printf(“%f\n”,DDA3[i]);

fprintf(fp3,“%f\n”,DDA3[i]);

//}

//else

//{

//printf(“%f\n”,DDA3[i]);

//fprintf(fp3,“%f\n”,DDA3[i]);

//}

}

for(i=0;i《aad3;i++)

{//if(i==0)

printf(“%s\n”,“AAD3”);

printf(“%f\n”,AAD3[i]);

fprintf(fp4,“%f\n”,AAD3[i]);

//}

//else

//{

//printf(“%f\n”,AAD3[i]);

//fprintf(fp4,“%f\n”,AAD3[i]);

//}

}

for(i=0;i《dad3;i++)

{//if(i==0)

printf(“%s\n”,“DAD3”);

printf(“%f\n”,DAD3[i]);

fprintf(fp5,“%f\n”,DAD3[i]);

//}

//else

//{

//printf(“%f\n”,DAD3[i]);

//fprintf(fp5,“%f\n”,DAD3[i]);

//}

}

for(i=0;i《add3;i++)

{//if(i==0)

printf(“%s\n”,“ADD3”);

printf(“%f\n”,ADD3[i]);

fprintf(fp6,“%f\n”,ADD3[i]);

//}

//else

//{

//printf(“%f\n”,ADD3[i]);

//fprintf(fp6,“%f\n”,ADD3[i]);

//}

//

}

for(i=0;i《ddd3;i++)

{//if(i==0)

printf(“%s\n”,“DDD3”);

printf(“%f\n”,DDD3[i]);

fprintf(fp7,“%f\n”,DDD3[i]);

//}

//else

//{

//printf(“%f\n”,DDD3[i]);

//fprintf(fp7,“%f\n”,DDD3[i]);

//}

}

//关闭文件

fclose(fp);

fclose(fp0);

fclose(fp1);

fclose(fp2);

fclose(fp3);

fclose(fp4);

fclose(fp5);

fclose(fp6);

fclose(fp7);

}

/*runthisprogramusingtheconsolepauseroraddyourowngetch,system(“pause”)orinputloop*/

/*尝试不对

doubleA1[(LENGTH+DB_LENGTH)/2];

doubleD1[(LENGTH+DB_LENGTH)/2];

inta1=(LENGTH+DB_LENGTH)/2;

intd1=(LENGTH+DB_LENGTH)/2;

doubleAA2[(LENGTH+DB_LENGTH)/4];

doubleDA2[(LENGTH+DB_LENGTH)/4];

inta2=(LENGTH+DB_LENGTH)/4;

intd2=(LENGTH+DB_LENGTH)/4;

doubleAAA3[(LENGTH+DB_LENGTH)/8];

doubleDAA3[(LENGTH+DB_LENGTH)/8];

inta3=(LENGTH+DB_LENGTH)/8;

intd3=(LENGTH+DB_LENGTH)/8;

doubleAAAA4[(LENGTH+DB_LENGTH)/16];

doubleDAAA4[(LENGTH+DB_LENGTH)/16];

Covlution(input,h,temp,n,m);

for(i=0;i《n+m-2;i+=2)

{

A1[i/2]=temp[i];//一层分解的低频部分,进行了2抽值,即尺度空间,低频概貌部分

}

Covlution(input,g,temp,n,m);

for(i=0;i《n+m-2;i+=2)

{

D1[i/2]=temp[i];//一层分解的高频部分,已经进行了2抽取,即高频细节部分

}

Covlution(A1,h,temp,a1,m);

for(i=0;i《a1+m-2;i+=2)

{

AA2[i/2]=temp[i];//一层分解的低频部分,进行了2抽值,即尺度空间,低频概貌部分

}

Covlution(A1,g,temp,a1,m);

for(i=0;i《a1+m-2;i+=2)

{

DA2[i/2]=temp[i];//一层分解的低频部分,进行了2抽值,即尺度空间,低频概貌部分

}

Covlution(AA2,h,temp,a2,m);

for(i=0;i《a2+m-2;i+=2)

{

AAA3[i/2]=temp[i];//一层分解的低频部分,进行了2抽值,即尺度空间,低频概貌部分

}

Covlution(AA2,g,temp,a2,m);

for(i=0;i《a2+m-2;i+=2)

{

DAA3[i/2]=temp[i];//一层分解的低频部分,进行了2抽值,即尺度空间,低频概貌部分

}

Covlution(AAA3,h,temp,a3,m);

for(i=0;i《a3+m-2;i+=2)

{

AAAA4[i/2]=temp[i];//一层分解的低频部分,进行了2抽值,即尺度空间,低频概貌部分

}

Covlution(AAA3,g,temp,a3,m);

for(i=0;i《a3+m-2;i+=2)

{

DAAA4[i/2]=temp[i];//一层分解的低频部分,进行了2抽值,即尺度空间,低频概貌部分

}

for(i=0;i《4116;i++)

{

if(i《=259)

output[i]=AAAA4[i];

elseif(i《=519)

output[i]=DAAA4[i-260];

elseif(i《=1035)

output[i]=DAA3[i-520];

elseif(i《=2064)

output[i]=DA2[i-1036];

else

output[i]=DA2[i-2065];

}

*/

打开APP阅读更多精彩内容
声明:本文内容及配图由入驻作者撰写或者入驻合作网站授权转载。文章观点仅代表作者本人,不代表电子发烧友网立场。文章及其配图仅供工程师学习之用,如有内容侵权或者其他违规问题,请联系本站处理。 举报投诉

全部0条评论

快来发表一下你的评论吧 !

×
20
完善资料,
赚取积分