#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;
}
转载于:https://my.oschina.net/u/582827/blog/1492209