小波包3層分解C語言源程序

ID:271983 · 發表于 2018-1-5 15:53

此程序是用于信號處理分析，突出奇異值的前段處理，對信號進行小波包分解，用C語言實現的，這個程序相比MATLAB程序，沒有進行邊緣的處理，直接根據卷積原理進行編寫而成的，如有對信號邊緣處理的要求，可以自行改進。程序在DEV 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];
}
*/

復制代碼

帳號		自動登錄	找回密碼
密碼			立即注冊

久久久久久久999_99精品久久精品一区二区爱城_成人欧美一区二区三区在线播放_国产精品日本一区二区不卡视频_国产午夜视频_欧美精品在线观看免费

小波包3層分解C語言源程序

評分