## 任意阶daubechies小波函数的构造 转

vga

``````#include <iostream.h>
#include <math.h>
#include <fstream.h>

#define Pi 3.14159265358968
#define N 512
void fft(double *x_real,double *x_image);
void ifft(double *x_real,double *x_image);
void complex_multip(double*,double*,double*,double*);
double* conv(double*,double*,int,int);
int exp2(int n);

void main()
{
ofstream file("E:/Debug/Mz.txt");

int p ;
cout<<"请输入所构造小波的消失矩：";
cin>>p;
cout<<endl;
double N_real[N],N_image[N],h[N];
double z_real[N],z_image[N],ones[N],Mz[N],Mz_image[N],temp1_real[N],temp1_image[N],temp2[N];
double G_real[N],G_image[N];

for(int k = 0; k < N ; k++)
{
z_real[k] = cos(2 * Pi * k / N);
z_image[k] = sin(2 * Pi * k / N);
ones[k] = 1.0;
Mz[k] = 0.0;
Mz_image[k] = 0.0;
temp1_real[k] = (1 + z_real[k])/2;
temp1_image[k] = -z_image[k]/2;
temp2[k] = (1 - z_real[k])/2;
N_real[k] = 1.0;
N_image[k] = 0.0;
}

for(int n = 0; n < p; n++)
{
for(k = 0 ; k < N ; k++)
{
Mz[k] += ones[k];
ones[k] = ones[k] * temp2[k] * (p + n) / (n + 1);
}
}

for(k = 0 ; k < N ; k++)
{
Mz[k] = 4 * Mz[k];
Mz[k] = log(Mz[k]);
}

ifft(Mz,Mz_image);

for(k = (N / 2) ; k < N; k++)
{
Mz[k] = 0.0;
Mz_image[k] = 0.0;
}
Mz[0] /= 2;
Mz_image[0] /= 2;

fft(Mz,Mz_image);

for(k = 0 ; k < N ; k++)
{
G_real[k] = exp(Mz[k]) * cos(Mz_image[k]);
G_image[k] = exp(Mz[k]) * sin(Mz_image[k]);
}

/*************************
计算(1+1/z)/2 的p次方

temp1 = (1+1/z)/2
N = ((1 + 1/z)/2)^p
***************************/
for(n = 0 ; n < p ; n++)
complex_multip(N_real,N_image,temp1_real,temp1_image);

/*G = N * G*/
complex_multip(G_real,G_image,N_real,N_image);

ifft(G_real,G_image);

for(k = 0 ; k < 2 * p ; k++)
h[k] = G_real[k]/sqrt(2);
for(k = 0 ; k < 2*p ; k++)
cout<<"h["<<k<<"]: "<<h[k]<<endl;
cout<<endl;
/      file<<h[k]<<' ';

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

以上由消失矩得到低通滤波器的系数

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

double FI[N],PSI[N];
for(int i = 0 ; i < 2 * p ; i++)
FI[i] = h[i];
int x_length = 2 * p ;
int y_length;

/*计算PHI(t)的系数*/
n = 1;
while(n<=2)
{
y_length = 2*p + (2*p-1) * (exp2(n) - 1);
double *y = new double[y_length];
double *x = new double[x_length+y_length-1];
int j = 0 ;
for(i = 0 ; i < y_length ; i++)
if(i % exp2(n) == 0)
y[i] = h[j++];
else
y[i] = 0;

x = conv(FI,y,x_length,y_length);
for(i = 0 ; i < x_length+y_length-1; i++)
FI[i] = x[i];
x_length = x_length + y_length - 1;

n++;
}

cout<<"下面是PHI(t)的系数："<<endl;
for(i = 0 ; i < x_length; i++)
cout<<FI[i]<<endl;
cout<<endl;

/* FI中存放PHI(t)的系数*/

/*计算PSI(t)的系数*/
double g[N];
for(i = 0 ; i < 2 * p ; i++)
{
if(i%2 == 0)
g[i] = h[2*p-i-1];
else
g[i] = -h[2*p-i-1];
PSI[i] = g[i];
}

x_length = 2 * p;
n = 1;
while(n<=2)
{
y_length = 2*p + (2*p-1) * (exp2(n) - 1);
double *y = new double[y_length];
double *x = new double[x_length+y_length-1];
int j = 0 ;
for(i = 0 ; i < y_length ; i++)
if(i % exp2(n) == 0)
y[i] = g[j++];
else
y[i] = 0;

x = conv(PSI,y,x_length,y_length);
for(i = 0 ; i < x_length+y_length-1; i++)
PSI[i] = x[i];
x_length = x_length + y_length - 1;

n++;
}

cout<<"下面是PSI(t)的系数："<<endl;
for(i = 0 ; i < x_length; i++)
cout<<PSI[i]<<endl;
}

void fft(double* x_real,double* x_image)
{
double *w_real,*w_image,t_real,t_image;
int i,j,k,n,n1;
int m;
int M = (int)(log(N)/log(2)); //分解的级数
double temp_real,temp_image;
w_real = new double[M];
w_image = new double[M];

n = 1;
for(m = 1 ; m < M ; m++)
{
n <<= 1;
w_real[m] = cos(Pi/n);
w_image[m] = -sin(Pi/n);
}

/*对数据预处理，交换位置*/
for(i = 0 ,j = 0; i < N-1 ; i++)
{
if(i < j)
{
temp_real = x_real[j];
x_real[j] = x_real[i];
x_real[i] = temp_real;

temp_image = x_image[j];
x_image[j] = x_image[i];
x_image[i] = temp_image;
}
k = N / 2;
while(k < (j + 1))
{
j -= k;
k >>= 1;
}
j += k;
}
/* m=0级，旋转因子为1 */
for(i = 0 ; i < N ; i+=2)
{
temp_real = x_real[i] - x_real[i+1];
temp_image = x_image[i] - x_image[i+1];
x_real[i] = x_real[i] + x_real[i+1];
x_image[i] = x_image[i] + x_image[i+1];
x_real[i+1] = temp_real;
x_image[i+1] = temp_image;
}
/*m=1，2，…，M级，旋转因子为W＝cos(Pi/2^m)-sin(Pi/2^m)*/
n = 2;
for(m = 1 ; m < M; m++)
{
n <<= 1;
n1 = n>>1;

t_real = 1.0;
t_image = 0.0;
for(j = 0 ; j < n1 ; j++)
{
for(i = j ; i < N ; i+=n)
{
k = i + n1;
temp_real = x_real[k] * t_real - x_image[k] * t_image;
temp_image = x_image[k] * t_real + x_real[k] * t_image;

x_real[k] = x_real[i] - temp_real;
x_image[k] = x_image[i] - temp_image;
x_real[i] = x_real[i] + temp_real;
x_image[i] = x_image[i] + temp_image;
}
temp_real = t_real;
t_real = t_real * w_real[m] - t_image * w_image[m];
t_image = temp_real * w_image[m] + t_image * w_real[m];
}
}
}

void ifft(double *x_real,double *x_image)
{
for(int i = 0 ; i < N; i++)
x_image[i] = -x_image[i];
fft(x_real,x_image);
for(i = 0 ; i < N; i++)
{
//      x_real[i] = -x_real[i];
x_image[i] = -x_image[i];
x_real[i] /= N;
x_image[i] /= N;
}
}

void complex_multip(double *x_real,double *x_image,double *y_real,double *y_image)
{
double t;
for(int k = 0 ; k < N; k++)
{
t = x_real[k] * y_real[k] - x_image[k] * y_image[k];
x_image[k] = x_real[k] * y_image[k] + x_image[k] * y_real[k];
x_real[k] = t;
}
}

double* conv(double *x,double *y,int x_length,int y_length)
{
//  double x_image[N],y_image[N];
double *temp = new double[x_length+y_length-1];
int i,j;
for(i = 0 ; i < x_length+y_length-1; i++)
temp[i] = 0.0;
for(i = 0 ; i < x_length + y_length -1; i++)
for(j = 0 ; j <= i; j++)
{
if( j >= x_length)
x[j] = 0;
if( i-j >= y_length)
y[i-j] = 0;
temp[i] += x[j] * y[i-j];
}

return temp;
}

int exp2(int n)
{
int result = 1;
for(int i = 0 ; i < n ; i++)
result *= 2;
return result;
}

``````

### 评论(7)

//读取输入信号
FILE *fp;
fp=fopen("data.txt","r");
if(fp==NULL) //如果读取失败
{
printf("错误！找不到要读取的文件/"data.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, nStep);

//一维小波变换后的结果写入txt文件
fp=fopen("test.txt","w");

//打印一维小波变换后的结果
for(i = 0; i < n/pow(2,nStep-1); i++)///pow(2,nStep-1)
{
printf("%f/n", data_output[i]);
fprintf(fp,"%f/n", data_output[i]);
}

//关闭文件
fclose(fp);
}

for(i = 0; i < nStep; i++)
{
Covlution(output, h, g, temp, n, m);
Sort(temp, output, n);
n = n/2;
}
}

void main()
{

double data[LENGTH];//输入信号
double temp[LENGTH];//中间结果
double data_output[LENGTH];//一维小波变换后的结果
int n = 0;//输入信号长度
int m = 6;//Daubechies正交小波基长度
int nStep = 6;//分解级数
int i = 0;
char s;//从txt文件中读取一行数据

static double h[] = {.332670552950, .806891509311, .459877502118,
-.135011020010, -.085441273882, .035226291882};
static double g[] = {.035226291882, .085441273882, -.135011020010,
-.459877502118, .806891509311, -.332670552950};

//读取输入信号

/*****************************************************************
*  排序函数
*
*  将卷积后的结果进行排序，使尺度系数和小波系数分开
*****************************************************************/
void Sort(double data[], double sort[], int n)
{
for(int i = 0; i < n; i+=2)
{
sort[i/2] = data[i];
}

for(i = 1; i < n; i+=2)
{
sort[n/2+i/2] = data[i];
}

}

/*****************************************************************
*　一维小波变换函数
*
*　说明: 一维小波变换,可进行多次分解
*
*　输入参数: input[],输入信号; output[],小波变换结果，包括尺度系数
*　和小波系数两部分; temp[],存放中间结果;h[],Daubechies小波基低通滤
*　波器系数;g[],Daubechies小波基高通滤波器系数;n,输入信号长度; m,
*　Daubechies小波基紧支集长度; nStep,小波变换分解次数
*
*　李承宇, lichengyu2345@126.com
*
* 2010-08-22
*****************************************************************/
void DWT1D(double input[], double output[], double temp[], double h[],
double g[], int n, int m, int nStep)
{
int i = 0;

for(i = 0; i < n; i++)
{
output[i] = input[i];
}

for(i = 0; i < nStep; i++)

//最后m/2-1行
//  i = ( (n - m) + m/2 + 1 ) ;//*********
for(j = 1; j <= m/2; j+=2, i+=2)
{
for(k = 0; k < j; k++)
{
cov[i] += data[k] * g[m-j-k];//k针对data[k]
}

for(k = 0; k < m-j; k++)
{
cov[i] += g[k] * data[n-(m-j)+k];//k针对core[k]
}
}
}

//****************************************************
//偶数行用g[]进行卷积
//****************************************************
//前m/2+1行
i = 1;
for(j = 0; j < m/2; j+=2, i+=2)
{
for(k = m/2-j; k < m; k++ )
{
cov[i] += data[k-(m/2-j)] * g[k];//k针对core[k]
}

for(k = n-m/2+j; k < n; k++ )
{
cov[i] += data[k] * g[k-(n-m/2+j)];//k针对data[k]
}
}

//中间的n-m行
for( ; i <= (n-m)+m/2; i+=2)
{
for( j = 0; j < m; j++)
{
cov[i] += data[i-m/2+j] * g[j];
}
}

//****************************************************
//奇数行用h[]进行卷积
//****************************************************
//前m/2+1行
i = 0;
for(j = 0; j < m/2; j+=2, i+=2)
{
for(k = m/2-j; k < m; k++ )
{
cov[i] += data[k-(m/2-j)] * h[k];//k针对core[k]
}

for(k = n-m/2+j; k < n; k++ )
{
cov[i] += data[k] * h[k-(n-m/2+j)];//k针对data[k]
}
}

//中间的n-m行
for( ; i <= (n-m)+m/2; i+=2)
{
for( j = 0; j < m; j++)
{
cov[i] += data[i-m/2+j] * h[j];
}
}

//最后m/2-1行
//  i = ( (n - m) + m/2 + 1 )/2*2;//**********
for(j = 1; j <= m/2; j+=2, i+=2)
{
for(k = 0; k < j; k++)
{
cov[i] += data[k] * h[m-j-k];//k针对data[k]
}

for(k = 0; k < m-j; k++)
{
cov[i] += h[k] * data[n-(m-j)+k];//k针对core[k]
}
}

1、题目：一维小波变换，可多次分解
2、原理：卷积核变为Daubechies正交小波基h[]和g[]的交替形式。增加了多次分解的功能。
3、代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define LENGTH 4096//信号长度
/*****************************************************************
*　一维卷积函数
*
*　说明: 循环卷积,卷积结果的长度与输入信号的长度相同
*
*　输入参数: data[],输入信号; h[],Daubechies小波基低通滤波器系数;
*　       g[],Daubechies小波基高通滤波器系数; cov[],卷积结果;
*       n,输入信号长度; m,卷积核长度.
*
*　李承宇, lichengyu2345@126.com
*
* 2010-08-22
*****************************************************************/
void Covlution(double data[], double h[], double g[], 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;
}

SLS机器学习介绍（05）：时间序列预测

00系列文章目录 0.1 算法原理目录 SLS机器学习介绍（01）：时序统计建模 SLS机器学习介绍（02）：时序聚类建模 SLS机器学习介绍（03）：时序异常检测建模 SLS机器学习介绍（04）：规则模式挖...

01/03
0
0

we_are_family678
2018/03/21
0
0
R语言数据挖掘实战系列（4）

R语言数据挖掘实战系列（4）——数据预处理 数据预处理一方面是要提高数据的质量，另一方面是要让数据更好地适应特定的挖掘技术或工具。数据预处理的主要内容包括数据清洗、数据集成、数据变...

PXZ6603
2017/07/12
0
0
4-1 Mallet小波分解重构程序

Mallet小波在小波届的地位类似fft在傅立叶变化中的地位，在分解过程中先滤波后抽取，重构过程中先插值后滤波，可以操作正交小波变换和双正交小波变换。 本文中的程序是对构造的信号进行高低通...

2014/04/14
2.9K
1

aresmiki
2017/02/17
0
0

var params = {};testParan.split(',').forEach(item =>{ var tmpArr = item.split('='); Vue.set(params, tmpArr[0].trim(), tmpArr[1].trim());});Vue.set(params, 'sql', sql);......

23分钟前
4
0

分立器件被广泛应用到消费电子、计算机及外设、网络通信，汽车电子、led显示屏等领域。 　　半导体产业中有两大分支：集成电路和分立器件。 　　集成电路 　　集成电路（integrated circ...

33分钟前
5
0
kibana rpm安装

36分钟前
3
0

37分钟前
4
0

SuspendThread()：挂起线程 If the function succeeds, the return value is the thread's previous suspend count; otherwise, it is (DWORD) -1. ResumeThread()：启动线程 If the functio......

rainbowcode
37分钟前
5
0