小波包3层分解C语言源程序 -

#include <stdio.h>
#include <stdlib.h>
#include <stdio.h>
#include <stdlib.h>
#define LENGTH 4096//信号长度
#define DB_LENGTH 8 //Daubechies小波基紧支集长度
/******************************************************************
*　一维卷积函数
*
*　说明: 循环卷积,卷积结果的长度与输入信号的长度相同
*
*　输入参数: data[],输入信号; core[],卷积核; cov[],卷积结果;
* n,输入信号长度; m,卷积核长度.
*
*　李承宇, lichengyu2345@126.com
*
* 2010-08-18
******************************************************************/
/*void Covlution(double data[], double core[], double cov[], int n, int m)
{
int i = 0;
int j = 0;
int k = 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]
}
}
}
*/
//定义一个线性卷积
void Covlution(double data[], double core[], double cov[], int n, int m)
{
int i = 0;
int j = 0;
int t = 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];
}
else if(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小波基紧支集长度.
*
*　李承宇, lichengyu2345@126.com
*
* 2010-08-19
******************************************************************/
void DWT1D_1(double input[], double output0[], double output1[],double temp[], double h[],
double g[], int n, int m)
{
// double temp[LENGTH] = {0};//?????????????
int i = 0;
/*
//尺度系数和小波系数放在一起
Covlution(input, h, temp, n, m);
for(i = 0; i < n; i += 2)
{
output[i] = te mp[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抽取，即高频细节部分
}
}
void DWT1D_2(double input[], double AA2[], double DA2[],double AD2[],double DD2[],double temp0[],double temp1[],double temp[], double h[],
double g[], int n, int m)
{
DWT1D_1(input, temp0, temp1,temp, h, g, n, m);
int a1=(m+n)/2;
DWT1D_1(temp0, AA2, DA2,temp, h, g, a1, m);
int d1=(n+m-4)/2;
DWT1D_1(temp1, AD2, DD2,temp, h, g, d1, m);
}
void DWT1D_3(double input[], double AAA3[], double DAA3[],double ADA3[],double DDA3[],double AAD3[], double DAD3[],double ADD3[],double DDD3[],
double temp0[],double temp1[],double temp2[],double temp3[],double temp00[], double temp11[], double temp[], double h[], double g[], int n, int m)
{
DWT1D_2(input, temp0, temp1,temp2, temp3,temp00, temp11,temp, h, g, n, m);
int aa2=(m+n)/4;
DWT1D_1(temp0, AAA3, DAA3,temp, h, g, aa2, m);
int da2=(n+3*m-8)/4;
DWT1D_1(temp1, ADA3, DDA3,temp, h, g, da2, m);
int ad2=(n+3*m-4)/4;
DWT1D_1(temp2, AAD3, DAD3,temp, h, g, ad2, m);
int dd2=(n+3*m-12)/4;
DWT1D_1(temp3, ADD3, DDD3,temp, h, g, dd2, m);
}
void main()
{
double data[LENGTH];//输入信号
double temp0[(LENGTH+DB_LENGTH)/4]; //先定义了一个中间结果
double temp1[(LENGTH+DB_LENGTH*3-8)/4];
double temp2[(LENGTH+DB_LENGTH*3-4)/4];
double temp3[(LENGTH+DB_LENGTH*3-12)/4];
double temp00[(LENGTH+DB_LENGTH)/2];
double temp11[(LENGTH+DB_LENGTH-4)/2];
double temp[LENGTH+DB_LENGTH-1];
double AAA3[(LENGTH+DB_LENGTH*5)/8];//一维小波变换后的结果
double DAA3[(LENGTH+DB_LENGTH*5-16)/8];
double ADA3[(LENGTH+DB_LENGTH*7-8)/8];
double DDA3[(LENGTH+DB_LENGTH*7-24)/8];
double AAD3[(LENGTH+DB_LENGTH*7-4)/8];
double DAD3[(LENGTH+DB_LENGTH*7-20)/8];
double ADD3[(LENGTH+DB_LENGTH*7-12)/8];
double DDD3[(LENGTH+DB_LENGTH*7-28)/8];
int aaa3=(LENGTH+DB_LENGTH*5)/8;//一维小波变换后的结果数组的长度
int daa3=(LENGTH+DB_LENGTH*5-16)/8;
int ada3=(LENGTH+DB_LENGTH*7-8)/8;
int dda3=(LENGTH+DB_LENGTH*7-24)/8;
int aad3=(LENGTH+DB_LENGTH*7-4)/8;
int dad3=(LENGTH+DB_LENGTH*7-20)/8;
int add3=(LENGTH+DB_LENGTH*7-12)/8;
int ddd3=(LENGTH+DB_LENGTH*7-28)/8;
int n = 0;//输入信号长度
int m = 8;//Daubechies正交小波基长度
int i = 0;
char s[32];//从txt文件中读取一行数据
/*//DB3
static double h[] = {.332670552950, .806891509311, .459877502118, -.135011020010,
-.085441273882, .035226291882};
static double g[] = {.035226291882, .085441273882, -.135011020010, -.459877502118,
.806891509311, -.332670552950};
*/
//DB4
static double h[] = {0.2303778133088964, 0.7148465705529154, 0.6308807679398587,-0.0279837694168599, -0.1870348117190931, 0.0308413818355607, 0.0328830116668852, -0.0105974017850690};//h[],Daubechies小波基低通滤波器系数;
static double g[] = {-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);
}
/* run this program using the console pauser or add your own getch, system("pause") or input loop */
/* 尝试不对
double A1[(LENGTH+DB_LENGTH)/2];
double D1[(LENGTH+DB_LENGTH)/2];
int a1=(LENGTH+DB_LENGTH)/2;
int d1=(LENGTH+DB_LENGTH)/2;
double AA2[(LENGTH+DB_LENGTH)/4];
double DA2[(LENGTH+DB_LENGTH)/4];
int a2=(LENGTH+DB_LENGTH)/4;
int d2=(LENGTH+DB_LENGTH)/4;
double AAA3[(LENGTH+DB_LENGTH)/8];
double DAA3[(LENGTH+DB_LENGTH)/8];
int a3=(LENGTH+DB_LENGTH)/8;
int d3=(LENGTH+DB_LENGTH)/8;
double AAAA4[(LENGTH+DB_LENGTH)/16];
double DAAA4[(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];
else if(i<=519)
output[i] = DAAA4[i-260];
else if(i<=1035)
output[i] = DAA3[i-520];
else if(i<=2064)
output[i] = DA2[i-1036];
else
output[i] = DA2[i-2065];
}
*/

复制代码