文档章节

任意阶daubechies小波函数的构造

vga
 vga
发布于 2017/07/30 05:50
字数 1061
阅读 56
收藏 0
#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;
}

本文转载自:来源不详

vga

vga

粉丝 23
博文 366
码字总数 26645
作品 0
佳木斯
私信 提问
加载中

评论(7)

vga
vga 博主
//读取输入信号
  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);
}
vga
vga 博主
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 s32;//从txt文件中读取一行数据

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

  //读取输入信号
vga
vga 博主

/*****************************************************************
*  排序函数
*
*  将卷积后的结果进行排序,使尺度系数和小波系数分开
*****************************************************************/
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++)
  
vga
vga 博主
//最后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]
    }
  }
}
vga
vga 博主
//****************************************************
  //偶数行用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];
    }
  }
vga
vga 博主
//****************************************************
  //奇数行用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]
    }
  }
vga
vga 博主
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
小波变换轻松入门(我的理解说明)

第一节 一个很简单的例子 还谈不上正式入门 但他具备了部分的思想。 [x0,x1,x2,x3]=[90,70,100,70] 为达到压缩 我们可取 (x0+x1)/2  (x0-x1)/2 来代表 x0,x1 这样 [90,70] 可表示为 [80,10...

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
希尔伯特变换和瞬时频率问题--连载(二)

写在开始的一段话: PS:OK,上一期关于希尔伯特变换的文章发出后,有知友在评论区说“看到最后……居然这……”,哈哈,其实我也挺愧疚大家的,明明一篇知识分享的文章,却写到结尾都没进入...

aresmiki
2017/02/17
0
0

没有更多内容

加载失败,请刷新页面

加载更多

将key=value转成对象形式

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安装

https://www.elastic.co/guide/en/kibana/6.2/rpm.html 下载对应的版本wget https://artifacts.elastic.co/downloads/kibana/kibana-6.2.4-x86_64.rpm 安装 rpm -ivh kibana-6.2.4-x86_64......

看的最远的地方
36分钟前
3
0
高防CDN相比较于高防服务器,为何更加稳定?

对于DDoS攻击,那些已经做过网站、平台的人应该知道,DDoS攻击是非常可怕的,因为这种攻击本质上不能防御,或者DDoS攻击只能被减轻,不能完全消除。DDoS,意思是“分布式拒绝服务”。它是一种...

云漫网络Ruan
37分钟前
4
0
线程SuspendThread() ResumeThread()的使用

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

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部