任意阶daubechies小波函数的构造

最新推荐文章于 2024-01-02 12:02:45 发布
最新推荐文章于 2024-01-02 12:02:45 发布
阅读量1.2k 收藏 4
点赞数
文章标签： python
原文链接：https://my.oschina.net/u/582827/blog/1492209
版权

2019独角兽企业重金招聘Python工程师标准>>> hot3.png

#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

weixin_34352005
关注
  • 0
    点赞
  • 踩
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
Daubechies小波介绍
weixin_34168880的博客
08-29 9224
Daubechies小波是正交、连续且紧支撑的。 正交条件下，$H(\omega)$必须满足下式： $|H(\omega)|^2+|H(\omega + \pi)|^2 =1$ 连续紧支撑条件下，$H(\omega)$必须满足下式： $H(\omega)=\left(  \frac{1+e^{-i\omega}}{2}  \right) ^N S(\omega)$ $S(\omega)=...
matlab中daubechie小波,[转载]小波滤波器
weixin_32538033的博客
04-06 1628
滤波器组完美重构与小波快速算法。前面的分析可以知道Vj相当于在j分辨率的逼近，Vj-1相当于j-1分辨率的逼近，这样Wj-1相当于两个分辨率逼近的差。在高分辨率下，我们可以用f在(2^j*t)的采样值来代替向Vj空间的投影，但是这是需要说明的，否则成为“小波的罪恶”，本来在Vj上的投影需要函数对Vj上的基{2^(j/2)*m(2^j*t–n)}投影，用采样值来代替是因为当j足够大的时候，如一般情况...
daubechies小波
12-05
daubechies小波的理论基础，介绍了daubechies小波的背景，构造结构，具有理论价值
使用matlab绘制N=2~7的Daubechies小波函数与尺度函数的波形图。
06-11
使用matlab绘制N=2~7的Daubechies小波函数与尺度函数的波形图。
小波的两个函数
ningyanggege的博客
07-17 2832
1、Haar小波         Haar，一般音译为“哈尔”。         Haar函数是小波分析中最早用到的一个具有紧支撑的正交小波函数，也是最简单的一个小波函数，它是支撑域在t∈[0,1]范围内的单个矩形波。         Haar小波在时域上是不连续的，所以作为基本小波性能不是特别好。   2、Daubechies(dbN)小波(紧支集正交小波)         Daub...
Daubechies 小波 Daubechies 小波 Daubechies 小波
03-29
对于较大的\(p\)值，尺度函数\(\varphi\)和小波函数\(\psi\)都表现出一致的Lipschitz正则性，即随着\(p\)的增加，它们将越来越接近于C^\infty级函数，这意味着它们不仅在任意阶导数上都是连续的，而且在任意点上都是...
区间B样条小波有限条法的构造研究
12-30
对于任意尺度j的m阶B样条尺度函数和相应的小波函数，可以通过特定的公式推导出来，进而用于有限条方法的构造。 7. 文章的资金支持信息：本研究得到高等学校博士学科点专项科研基金资助项目和上海教育委员会资助项目...
小波概念理解
leemboy的博客
10-30 6281
Daubechies小波 支撑长度         小波函数Ψ(t)、Ψ(ω)、尺度函数φ(t)和φ(ω)的支撑区间，是当时间或频率趋向于无穷大时，Ψ(t)、Ψ(ω)、φ(t)和φ(ω)从一个有限值收敛到0的长度。支撑长度越长，一般需要耗费更多的计算时间，且产生更多高幅值的小波系数。大部分应用选择支撑长度为5~9之间的小波，因为支撑长度太长会产生边界问题，支撑长度太短消失矩太低，不利于信号能...
小波变换中的多贝西小波（DB小波函数）概述
热门推荐
ASS-ASH的博客
07-21 1万+
英语DaubechiesWavelet），是以比利时女性物理暨数学家（IngridDaubechies）的名字命名之一种函数，当初英格丽·多贝西发现了一种具有阶层（hierarchy）性质的小波，便将此小波以她的名字命名。多贝西小波主要应用在离散型的，是最常使用到的小波转换，通常使用在、跟去除。一般而言的离散小波转换通常是以（orthogonalwavelet）为基底，而多贝西小波也是一种正交小波。由于它很容易经由的系数，称作（scalingfilter,也称为fatherwavelet）......
python小波阈值去噪
08-06
利用小波通用阈值对离散数据集进行滤波，程序中小波基和阈值函数可改
特征变换（4）Daubechies小波变换
Rana专栏
03-17 8908
Daubechies4小波变换可定义为下列变换矩阵 . 空白处表示零。c0、c1、c2和c3必须满足如下正交属性 c02 + c12 + c22 + c32 = 1 c2c0 + c3c1 = 0 c3 – c2 + c1 – c0 = 0 0c3 – 1c2 + 2c1 – 3c0 = 0 有唯一解 . 数组X的Daubechies4小波变换
用迭代法求Daubechies小波的尺度函数
zhaopeng
01-07 1088
  p = [0.683013,1.183013,0.316987,-0.183013]; N = 20; phi_before = 1; for i = 1:N;     phi_after = conv(p,phi_before);     n = length(phi_after);     phi_before(1:2:2*n) = phi_after;
python 用二维对称Daubechies小波对图像变换demo程序
leemboy的博客
11-30 3559
用Daubeches db2对称小波变化的lena图像 详细python代码如下，使用格式建代码说明部分 #this demo for 2D wavelets transform is developed by Xueshang Nai #Referrence Link:https://pywavelets.readthedocs.io/en/latest/regression/wp2d...
小波分析——2. 小波函数及其逆函数
老程的技术笔记
01-07 5446
由于小波变换作为一种可以将原始信号分解为不同频带的数学工具，而且在上一章中我们已经简要的介绍了它与傅里叶函数的异同，和一些其他特点，所以在这一章中我们再继续探讨一些关于小波在信号分解和合成方面的数学原理。Coiflets 小波是一类常用的小波核函数，具有较高的时域精度和较低的频域精度。所以，我们可以立即得到一个和上面相似的，但是更简单点的方法来描述这组公式，也就是小波的离散形式，即对于离散信号。Symlets 小波是一类常用的小波核函数，其中最常用的是。Haar小波是最简单的小波，其公式为。
python小波分解1
最新发布
weixin_44245188的博客
01-02 1092
安装包是PyWavelets，不是pywt，只是引入名字叫pywt，可以使用python -m pip install PyWavelets命令安装。
【信号去噪】基于小波变换实现脉搏信号去噪附Matlab代码
Matlab_dashi的博客
09-21 704
小波变换可用于脉搏信号去噪的目的。脉搏信号是一种包含大量噪声的生物信号，而小波变换可以通过将信号分解成不同频率的细节和近似部分来分析和处理信号。通过对细节部分进行阈值处理或滤波操作，可以去除噪声并保留有用的信号信息。在小波变换中，信号首先与一组特定的小波基函数进行卷积，然后通过缩放和平移来表示不同频率的分量。通过选择适当的小波基函数和设置合适的阈值，可以消除脉搏信号中的噪声。
小波分析的理解
weixin_30276935的博客
08-08 4279
　　小波变换是克服其他信号处理技术缺陷的一种分析信号的方法。 　　小波由一族小波基函数构成，它可以描述信号时间（空间）和频率（尺度）域的局部特性。采用小波分析最大优点是可对信号进行实施局部分析，可在任意的时间或空间域中分析信号。小波分析具有发现其他信号分析方法所不能识别的、隐藏于数据之中的表现结构特性的信息，而这些特性对机械故障和材料的损伤等识别是尤为重要的。如何选择小波基函数目前还没有一个理论...
python实现信号的抽样和内插
whujk的博客
09-11 2501
###matlab转python实现系列（二） 信号的抽样和内插 程序结构 #mermaid-svg-qOdTSPkGw8OtHn8i .label{font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family);fill:#333;color:#333}#mermaid-svg-qOdTSPkGw8OtHn8i .label text{fill:#333}#mermaid-svg-qOdTSPkGw8OtH
用matlab绘制三阶daubechies小波函数图像
04-19
好的，这是一个技术问题，我可以回答。以下是Matlab代码示例： ```matlab % 设置小波系数 h = [ (1+sqrt(3))/(4*sqrt(2)), (3+sqrt(3))/(4*sqrt(2)), (3-sqrt(3))/(4*sqrt(2)), (1-sqrt(3))/(4*sqrt(2))]; % 设置小波函数的区间 t = 0:0.001:1; % 计算小波函数值 w = zeros(size(t)); for i=0:length(h)-1 w = w + h(i+1).*phi((2*t-i)/2); end % 绘制小波函数图像 plot(t, w); title('Three-order Daubechies Wavelet'); xlabel('Time (s)'); ylabel('Amplitude'); ``` 其中，phi(x)为Haar函数，实现代码如下： ```matlab function y = phi(x) y = (x>=0) & (x<1); end ```
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值