离散小波变换 订阅
离散小波变换是对基本小波的尺度和平移进行离散化。在图像处理中,常采用二进小波作为小波变换函数,即使用2的整数次幂进行划分。 展开全文
离散小波变换是对基本小波的尺度和平移进行离散化。在图像处理中,常采用二进小波作为小波变换函数,即使用2的整数次幂进行划分。
信息
外文名
Discrete wavelettransform
中文名
离散小波变换
直方图定义
余弦变换是经典的谱分析工具,他考察的是整个时域过程的频域特征或整个频域过程的时域特征,因此对于平稳过程,他有很好的效果,但对于非平稳过程,他却有诸多不足。在JPEG中,离散余弦变换将图像压缩为8×8 的小块,然后依次放入文件中,这种算法靠丢弃频率信息实现压缩,因而图像的压缩率越高,频率信息被丢弃的越多。在极端情况下,JPEG图像只保留了反映图像外貌的基本信息,精细的图像细节都损失了。小波变换是现代谱分析工具,他既能考察局部时域过程的频域特征,又能考察局部频域过程的时域特征,因此即使对于非平稳过程,处理起来也得心应手。他能将图像变换为一系列小波系数,这些系数可以被高效压缩和存储,此外,小波的粗略边缘可以更好地表现图像,因为他消除了DCT压缩普遍具有的方块效应。
收起全文
精华内容
参与话题
问答
  • 离散小波变换DWT

    2019-04-08 10:00:38
    https://blog.csdn.net/u012507022/article/details/50979009)强推
    展开全文
  • 快速傅氏变换和离散小波变换附Matlab程序-小波变换 dwt计算方法 sy10.doc 长期以来,快速傅氏变换和离散小波变换在数字信号处理、石油勘探、地震预报、医学断层诊断、编码理论、量子物理及概率论等领域中都得到了...
  • 离散序列的Mallat算法分解公式如下: 其中,H(n)、G(n)分别表示所选取的小波函数对应的低通和高通滤波器的抽头系数序列。 从Mallat算法的分解原理可知,分解后的序列就是原序列与滤波器序列的卷积再进行隔点...

    1 Mallat算法

    离散序列的Mallat算法分解公式如下:

    其中,H(n)、G(n)分别表示所选取的小波函数对应的低通和高通滤波器的抽头系数序列。

    从Mallat算法的分解原理可知,分解后的序列就是原序列与滤波器序列的卷积再进行隔点抽取而来。

     

    离散序列的Mallat算法重构公式如下:

    其中,h(n)、g(n)分别表示所选取的小波函数对应的低通和高通滤波器的抽头系数序列。


    2 小波变换实现过程(C/C++)

    2.1       小波变换结果序列长度

          小波的Mallat算法分解后的序列长度由原序列长SoureLen和滤波器长FilterLen决定。从Mallat算法的分解原理可知,分解后的序列就是原序列与滤波器序列的卷积再进行隔点抽取而来。即分解抽取的结果长度为(SoureLen+FilterLen-1)/2。

    2.2       获取滤波器组

    对于一些通用的小波函数,简单起见,可以通过Matlab的wfilters(‘wavename’)获取4个滤波器;特殊的小波函数需要自行构造获得。

    下面以db1小波函数(Haar小波)为例,其变换与重构滤波器组的结果如下:

    //matlab输入获取命令
    >> [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('db1')
    
    //获取的结果
    Lo_D =
        0.7071    0.7071
    Hi_D =
       -0.7071    0.7071
    Lo_R =
        0.7071    0.7071
    Hi_R =
        0.7071   -0.7071
    
    2.3       信号边界延拓

            在Mallat算法中,假定输入序列是无限长的,而实际应用中输入的信号是有限的采样序列,这就会出现信号边界处理问题。对于边界信号的延拓一般有3种方法,即零延拓、对称延拓和周期延拓。

    3种延拓方法比较情况如下:

            对于正交小波变换来说,前两种延拓方法实现起来比较简单,但重建时会产生边界效应,而且分解的层数越多,产生的边界效应越显著。零延拓方法给人一种跳跃的感觉。至于对称性延拓,由于正交小波滤波器一般都是非对称性的(Haar小波基虽然是正交的,但它是非连续的),重建图象给人一种错位的感觉。相比较而言,只有最后一种延拓方式可以得到比较精确的重建结果,它不仅能保证分解与重建正确计算,而且恢复的质量也好。不过,周期性延拓方法虽然是常用的三种方法中比较好的方法,但会导致信号边缘的非连续性,从而会使得较高频率(子带)层的小波系数很大,即使信号本身相当平滑。从信号压缩的角度看,大的系数是希望避免的。信号的对称延拓可避免边缘的非连续性问题。然而,对称延拓只能和对称的小波滤波器一起适用。如果降低正交性要求,选择双正交小波变换,对称性延拓不失为一种好的方法。周期延拓可适用于任何小波变换,但可能导致输入序列边缘的不连续,使得高频系数较大。而对称延拓则避免了输入序列边界的不连续,是当前广泛采用的一种延拓方法。

    下式中给出了最常用的对称延拓表达式。



    当原序列长sLEN为偶数时延拓后的序列长为sLEN+2*(filterLEN),而原序列长为奇数时则需要在右端再延拓一个元素。注:在Matlab中默认使用了对称延拓。

    2.4       小波分解

    在db1小波函数下,离散序列的Mallat算法分解公式展开如下:


    其它的db小波,不再详述。小波分解C++源码如下。

    /**
     * @brief 小波变换之分解
     * @param sourceData 源数据
     * @param dataLen 源数据长度
     * @param db 过滤器类型
     * @param cA 分解后的近似部分序列-低频部分
     * @param cD 分解后的细节部分序列-高频部分
     * @return
     *         正常则返回分解后序列的数据长度,错误则返回-1
     */
    int Wavelet::Dwt(double *sourceData, int dataLen, Filter db, double *cA, double *cD)
    {
        if(dataLen < 2)
            return -1;
        if((NULL == sourceData)||(NULL == cA)||(NULL == cD))
            return -1;
    
        m_db = db;
        int filterLen = m_db.length;
        int n,k,p;
        int decLen = (dataLen+filterLen-1)/2;
        double tmp=0;
        cout<<"decLen="<<decLen<<endl;
    
        for(n=0; n<decLen; n++)
        {
            cA[n] = 0;
            cD[n] = 0;
            for(k=0; k<filterLen; k++)
            {
                p = 2*n-k;  // 感谢网友onetrain指正,此处由p = 2*n-k+1改为p = 2*n-k,解决小波重构后边界异常问题--2013.3.29 by hmm
    
                // 信号边沿对称延拓
                if((p<0)&&(p>=-filterLen+1))
                    tmp = sourceData[-p-1];
                else if((p>dataLen-1)&&(p<=dataLen+filterLen-2))
                    tmp = sourceData[2*dataLen-p-1];
                else if((p>=0)&&(p<dataLen-1))
                    tmp = sourceData[p];
                else
                    tmp = 0;
    
                // 分解后的近似部分序列-低频部分
                cA[n] += m_db.lowFilterDec[k]*tmp;
    
                // 分解后的细节部分序列-高频部分
                cD[n] += m_db.highFilterDec[k]*tmp;
            }
    
        }
    
        return decLen;
    }
    

    2.5      小波重构
    /**
     * @brief 小波变换之重构
     * @param cA 分解后的近似部分序列-低频部分
     * @param cD 分解后的细节部分序列-高频部分
     * @param cALength 输入数据长度
     * @param db 过滤器类型
     * @param recData 重构后输出的数据
     */
    void  Wavelet::Idwt(double *cA, double *cD,  int cALength, Filter db, double *recData)
    {
        if((NULL == cA)||(NULL == cD)||(NULL == recData))
            return;
    
        m_db = db;
        int filterLen = m_db.length;
    
        int n,k,p;
        int recLen = 2*cALength-filterLen+1;
        cout<<"recLen="<<recLen<<endl;
    
        for(n=0; n<recLen; n++)
        {
            recData[n] = 0;
            for(k=0; k<cALength; k++)
            {
                p = n-2*k+filterLen-1;
    
                // 信号重构
                if((p>=0)&&(p<filterLen))
                {
                    recData[n] += m_db.lowFilterRec[p]*cA[k] + m_db.highFilterRec[p]*cD[k];
                    //cout<<"recData["<<n<<"]="<<recData[n]<<endl;
                }
    
            }
        }
    }
    
    2.6      c++实现结果

    3 小波变换实现(MATLAB)

    使用matlab小波工具箱实现db4的情况如下。

    1、MatlabDB4.m文件内容。

    %加载txt数据示例
    s=importdata('data2.txt'); %load data2.txt;
    subplot(521);plot(s); %画出原始信号的波形图
    title('原始信号');
    
    [cA,cD]=dwt(s,'db4'); %采用db4小波并对信号进行一维离散小波分解。
    y=idwt(cA,cD,'db4'); %一维离散小波反变换
    subplot(522);
    plot(cA); %画出波形图
    title('MATLAB低频部分dwt-cA');
    
    subplot(523);
    plot(cD); %画出波形图
    title('MATLAB高频部分dwt-cD');
    
    subplot(524);
    plot(y); %画出波形图
    title('MATLAB重构idwt');
    

    2、波形图如下。

    4 小结

            在此,采用C++实现了dbN系列小波的单一维离散小波变换,通过对比db4小波变换发现(笔者也同样对比验证了db1-db3),C++实现效果和Matlab效果完全一样。通过这一过程,可以很好地帮助理解小波变换的Mallat算法原理,并将C++代码快速应用到工程实践,同时对MATLAB小波工具箱的实现细节也有更深入的理解。相关文献表明,C++代码实现的DWT算法比Matlab小波工具箱dwt方法效率更高。


    注:对小波变换见识还不深,有问题者,欢迎提出讨论交流,对源码细节感兴趣的tx,可以到如下链接下载。

    http://download.csdn.net/detail/v_hyx/5113111

    尺度
    展开全文
  • 在上回《小波学习之一》中,已经详细介绍了Mallat算法... MATLAB作为经典的数学工具,分析其小波变换dwt和idwt实现后发现真的很经典,学习参考价值很高。下面结合南京理工大学 谭彩铭的《解读matlab之小波库函数》及

            在上回《小波学习之一》中,已经详细介绍了Mallat算法C++实现,效果还可以,但也存在一些问题,比如,代码难于理解,同时出现了边界问题。在此,本文将重构代码,采用新的方法解决这些问题,同时也加深对小波变换的理解。

            MATLAB作为经典的数学工具,分析其小波变换dwt和idwt实现后发现真的很经典,学习参考价值很高。下面结合南京理工大学 谭彩铭的《解读matlab之小波库函数》及MATLAB小波工具包中m文件的情况,作一个小结,最后用C++函数进行实现,并且编译调试OK。
            一、MATLAB上dwt函数的工作过程
            假设x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)],计算y=dwt(x,’db2’),其计算过程主要由三个部分组成:
    1、边缘延拓,它主要由函数wextend完成。

    仔细分析子程序部分,函数wextend的用法为y=wextend('1D','sym',x,3);这样得到的y=[ x(3) x(2) x(1) x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(7) x(6) x(5)]
    2、卷积运算,它主要由函数conv2完成。
    仔细分析子程序部分,核心语句有z=conv2(y,Lo_D,'valid');这里设Lo_D=[h(1) h(2) h(3) h(4)]。

    这2步的实现过程示意图如下:



    3、最后就是下采样即隔点采样,其下采样是按照式a = z(2:2:length(z))进行的,高频低频部分均如此,项数为floor((7+4-1)/2)。

    最后的dwt低频系数结果是[z(2) z(4) z(6) z(8) z(10)],高频系数求解过程和低频系数一样,在此不再赘述。


            二、MATLAB上idwt函数的工作过程

    1、上采样即隔点插0,dyadup(x,0)。

    2、卷积运算,它也是最终由函数conv2完成。

    3、抽取结果,wkeep1(x,s,'c')。


    下面啥都不说show核心代码实现,欢迎讨论。

    /**
     * @brief 边缘延拓
     * @param typeId 延拓数据的类型,1D or 2D
     * @param modeId 延拓方式:对称、周期
     * @param in 输入数据
     * @param inLen 输入数据的长度
     * @param filterLen 小波基滤波器长度
     * @param out 返回结果数组
     * @return 返回结果数组长度
     */
    int SignalExtension(int typeId,
    		int modeId,   
    		double *in,   
    		int inLen,    
    		int filterLen, 
    		double out[])  
    {
        if((NULL == in)||(NULL == out))
            return -1;
        if(0 != typeId) // 目前只支持一种模型
        	return -1;
        //if(0 != modeId) // 目前只支持一种模型,信号对称拓延  'sym' or 'symh'  	Symmetric-padding (half-point): boundary value symmetric replication
        //	return -1;
        if( inLen < filterLen ) // inLen should lager than or equal extendLen, otherwise no extension
        	return -1;
    
        int i;
        int extendLen = filterLen - 1;
    
        if(0 == modeId) // 信号对称拓延
        {
            for(i=0; i<inLen; i++)
            {
            	out[extendLen+i] = in[i];
            }
            for(i=0; i<extendLen; i++)
            {
            	out[i]                     = out[2*extendLen - i - 1];       // 左边沿对称延拓
            	out[inLen + extendLen + i] = out[extendLen + inLen - i - 1]; // 右边沿对称延拓
            }
    
            return inLen + 2*extendLen;
        }
        else if(1 == modeId) // 信号周期拓延
        {
    		for( i = 0; i < extendLen; i++ )
    			out[i] = in[inLen-extendLen+i];
    		for ( i = 0; i < inLen; i++ )
    			out[extendLen+i] = in[i];
    
            return inLen + extendLen;
        }
    
    }

    /**
     * @brief 上采样  隔点插0
     * @param data 输入数据指针
     * @param n 输入数据长度
     * @param result 返回结果数组
     * @return 返回结果数组长度
     */
    int Upsampling(double* data, int n, double result[])
    {
    	int i;
    
    	for( i = 0; i < n; i++ )
    	{
    		result[2*i] = data[i];
    		result[2*i+1] = 0;
    	}
    
    	return( 2*n );
    }
    /**
     * @brief 下采样  隔点采样
     * @param data 输入数据指针
     * @param n 输入数据长度
     * @param result 返回结果数组
     * @return 返回结果数组长度
     */
    int Downsampling(double* data, int n, double result[])
    {
    	int i, m;
    
    	m = n/2;
    	for( i = 0; i < m; i++ )
    		result[i] = data[i*2 + 1];
    
    	return( m );
    }


    /**
     * @brief 卷积运算
     * @param shapeId 卷积结果处理方式
     * @param double *inSignal, int signalLen, // 输入信号及其长度
     * @param double *inFilter, int filterLen, // 输入滤波器及其长度
     * @param double outConv[], int *convLen)   // 输出卷积结果及其长度
     * @return
     */
    void Conv1(int shapeId,                  // 卷积结果处理方式
    		double *inSignal, int signalLen, // 输入信号及其长度
    		double *inFilter, int filterLen, // 输入滤波器及其长度
    		double outConv[], int *convLen)   // 输出卷积结果及其长度
    {
        if((NULL == inSignal)||(NULL == inFilter)||(NULL == outConv))
            return;
    
        int n,k,kmin,kmax,p;
        if(0 == shapeId)      // 对于MATLAB conv(...,'shape')  -----full
        {
        	*convLen = signalLen + filterLen - 1;
        	for (n = 0; n < *convLen; n++)
        	{
        		outConv[n] = 0;
    
        	    kmin = (n >= filterLen - 1) ? n - (filterLen - 1) : 0;
        	    kmax = (n < signalLen - 1) ? n : signalLen - 1;
    
        	    for (k = kmin; k <= kmax; k++)
        	    {
        	    	outConv[n] += inSignal[k] * inFilter[n - k];
        	    }
        	}
        }
        else if(1 == shapeId) // 对于MATLAB conv(...,'shape')  -----valid
        {
        	*convLen = signalLen - filterLen + 1;
        	for (n = filterLen - 1; n < signalLen; n++)
        	{
        		p = n - filterLen + 1;
        		outConv[p] = 0;
    
        	    kmin = (n >= filterLen - 1) ? n - (filterLen - 1) : 0;
        	    kmax = (n < signalLen - 1) ? n : signalLen - 1;
    
        	    for (k = kmin; k <= kmax; k++)
        	    {
        	    	outConv[p] += inSignal[k] * inFilter[n - k];
        	    }
        	}
        }
        else
        	return ;
    
    }

    /**
     * @brief 小波变换之分解
     * @param sourceData 源数据
     * @param dataLen 源数据长度
     * @param db 过滤器类型
     * @param cA 分解后的近似部分序列-低频部分
     * @param cD 分解后的细节部分序列-高频部分
     * @return 正常则返回分解后序列的数据长度,错误则返回-1
     */
    int Wavelet::Decomposition(double* sourceData, int dataLen, Filter db, double* cA, double* cD)
    {
        if(dataLen < 2)
            return -1;
        if((NULL == sourceData)||(NULL == cA)||(NULL == cD))
            return -1;
    
        m_db = db;
        int filterLen = m_db.length;
        int i, n;
        int decLen = (dataLen+filterLen-1)/2;
        int convLen = 0;
        double extendData[dataLen+2*filterLen-2];
        double convDataLow[dataLen+filterLen-1];
        double convDataHigh[dataLen+filterLen-1];
    
    /*
    MATLAB上dwt函数的工作过程
    假设x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)],计算y=dwt(x,’db2’)。
    其计算过程主要由两个部分组成:
    1:边缘延拓,它主要由函数wextend完成。
    2:卷积运算,它主要由函数conv2完成。
    先看第一部分,仔细分析子程序部分,函数wextend的用法为y=wextend('1D','sym',x,3);
    这样得到的y=[ x(3) x(2) x(1) x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(7) x(6) x(5)]
    在看第二部分,仔细分析子程序部分,核心语句有z=conv2(y,Lo_D,'valid');
    这里设Lo_D=[h(1) h(2) h(3) h(4)]。
    3:最后就是下采样,其下采样是按照式a = z(2:2:length(z))进行的,高频低频部分均如此,项数为floor((7+4-1)/2)。
     */
        // 1.边缘延拓
        SignalExtension(0, 0 , sourceData, dataLen, filterLen, extendData);
    
        // 2.卷积运算
        Conv1(1, extendData, dataLen+2*filterLen-2, db.lowFilterDec, filterLen, convDataLow, &convLen);
        Conv1(1, extendData, dataLen+2*filterLen-2, db.highFilterDec, filterLen, convDataHigh, &convLen);
    
        // 3.下采样
        Downsampling(convDataLow, dataLen + filterLen - 1, cA);
        Downsampling(convDataHigh, dataLen + filterLen - 1, cD);
    
        return decLen;
    }


    /**
     * @brief 小波变换之重构
     * @param cA 分解后的近似部分序列-低频部分
     * @param cD 分解后的细节部分序列-高频部分
     * @param cALength 输入数据长度
     * @param RecLength 输入重构后的原始数据长度
     * @param db 过滤器类型
     * @param recData 重构后输出的数据
     * @return 正常则返回重构数据长度,错误则返回-1
     */
    int Wavelet::Reconstruction(double *cA, double *cD, int cALength, int RecLength, Filter db, double* recData)
    {
        if((NULL == cA)||(NULL == cD)||(NULL == recData))
            return -1;
    
        m_db = db;
        int filterLen = m_db.length;
    
        int i,j;
        int n,k,p;
        int recLen = RecLength;
    
        int convLen = 0;
        double convDataLow[recLen+filterLen-1];
        double convDataHigh[recLen+filterLen-1];
    
        double cATemp[2*cALength];
        double cDTemp[2*cALength];
    
        memset(convDataLow, 0, (recLen+filterLen-1)*sizeof(double)); // 清0
        memset(convDataHigh, 0, (recLen+filterLen-1)*sizeof(double)); // 清0
        memset(cATemp, 0, 2*cALength*sizeof(double)); // 清0
        memset(cDTemp, 0, 2*cALength*sizeof(double)); // 清0
    
        // 1.隔点插0
        Upsampling(cA, cALength, cATemp);
        Upsampling(cD, cALength, cDTemp);
    
        // 2.卷积运算
        Conv1(0, cATemp, 2*cALength-1, db.lowFilterRec, filterLen ,convDataLow, &convLen);
        convLen = 0;
        Conv1(0, cDTemp, 2*cALength-1, db.highFilterRec, filterLen ,convDataHigh, &convLen);
    
        // 3.抽取结果及求和——实现类似MATLAB中的wkeep1(s,len,'c')的功能
        k = (convLen - recLen)/2;
        for(i=0; i<recLen; i++)
        {
        	recData[i] = convDataLow[i + k] + convDataHigh[i + k];
        }
    	
        return recLen;
    }



    尺度
    展开全文
  • 小波学习之二(单层一维离散小波变换DWT的Mallat算法C++实现优化) 在上回《小波学习之一》中,已经详细介绍了Mallat算法C++实现,效果还可以,但也存在一些问题,比如,代码难于理解,同时出现了边界问题。...

    小波学习之二(单层一维离散小波变换DWT的Mallat算法C++实现优化)

     

            在上回《小波学习之一》中,已经详细介绍了Mallat算法C++实现,效果还可以,但也存在一些问题,比如,代码难于理解,同时出现了边界问题。在此,本文将重构代码,采用新的方法解决这些问题,同时也加深对小波变换的理解。

            MATLAB作为经典的数学工具,分析其小波变换dwt和idwt实现后发现真的很经典,学习参考价值很高。下面结合南京理工大学 谭彩铭的《解读matlab之小波库函数》及MATLAB小波工具包中m文件的情况,作一个小结,最后用C++函数进行实现,并且编译调试OK。
            一、MATLAB上dwt函数的工作过程
            假设x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)],计算y=dwt(x,’db2’),其计算过程主要由三个部分组成:
    1、边缘延拓,它主要由函数wextend完成。

    仔细分析子程序部分,函数wextend的用法为y=wextend('1D','sym',x,3);这样得到的y=[ x(3) x(2) x(1) x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(7) x(6) x(5)]
    2、卷积运算,它主要由函数conv2完成。
    仔细分析子程序部分,核心语句有z=conv2(y,Lo_D,'valid');这里设Lo_D=[h(1) h(2) h(3) h(4)]。

    这2步的实现过程示意图如下:


    3、最后就是下采样即隔点采样,其下采样是按照式a = z(2:2:length(z))进行的,高频低频部分均如此,项数为floor((7+4-1)/2)。

    最后的dwt低频系数结果是[z(2) z(4) z(6) z(8) z(10)],高频系数求解过程和低频系数一样,在此不再赘述。

     

            二、MATLAB上idwt函数的工作过程

    1、上采样即隔点插0,dyadup(x,0)。

    2、卷积运算,它也是最终由函数conv2完成。

    3、抽取结果,wkeep1(x,s,'c')。

     

    下面啥都不说show核心代码实现,欢迎讨论。

     

    1. /** 
    2.  * @brief 边缘延拓 
    3.  * @param typeId 延拓数据的类型,1D or 2D 
    4.  * @param modeId 延拓方式:对称、周期 
    5.  * @param in 输入数据 
    6.  * @param inLen 输入数据的长度 
    7.  * @param filterLen 小波基滤波器长度 
    8.  * @param out 返回结果数组 
    9.  * @return 返回结果数组长度 
    10.  */  
    11. int SignalExtension(int typeId,  
    12.         int modeId,     
    13.         double *in,     
    14.         int inLen,      
    15.         int filterLen,   
    16.         double out[])    
    17. {  
    18.     if((NULL == in)||(NULL == out))  
    19.         return -1;  
    20.     if(0 != typeId) // 目前只支持一种模型  
    21.         return -1;  
    22.     //if(0 != modeId) // 目前只支持一种模型,信号对称拓延  'sym' or 'symh'      Symmetric-padding (half-point): boundary value symmetric replication  
    23.     //  return -1;  
    24.     if( inLen < filterLen ) // inLen should lager than or equal extendLen, otherwise no extension  
    25.         return -1;  
    26.   
    27.     int i;  
    28.     int extendLen = filterLen - 1;  
    29.   
    30.     if(0 == modeId) // 信号对称拓延  
    31.     {  
    32.         for(i=0; i<inLen; i++)  
    33.         {  
    34.             out[extendLen+i] = in[i];  
    35.         }  
    36.         for(i=0; i<extendLen; i++)  
    37.         {  
    38.             out[i]                     = out[2*extendLen - i - 1];       // 左边沿对称延拓  
    39.             out[inLen + extendLen + i] = out[extendLen + inLen - i - 1]; // 右边沿对称延拓  
    40.         }  
    41.   
    42.         return inLen + 2*extendLen;  
    43.     }  
    44.     else if(1 == modeId) // 信号周期拓延  
    45.     {  
    46.         for( i = 0; i < extendLen; i++ )  
    47.             out[i] = in[inLen-extendLen+i];  
    48.         for ( i = 0; i < inLen; i++ )  
    49.             out[extendLen+i] = in[i];  
    50.   
    51.         return inLen + extendLen;  
    52.     }  
    53.   
    54. }  
    /**
     * @brief 边缘延拓
     * @param typeId 延拓数据的类型,1D or 2D
     * @param modeId 延拓方式:对称、周期
     * @param in 输入数据
     * @param inLen 输入数据的长度
     * @param filterLen 小波基滤波器长度
     * @param out 返回结果数组
     * @return 返回结果数组长度
     */
    int SignalExtension(int typeId,
    		int modeId,   
    		double *in,   
    		int inLen,    
    		int filterLen, 
    		double out[])  
    {
        if((NULL == in)||(NULL == out))
            return -1;
        if(0 != typeId) // 目前只支持一种模型
        	return -1;
        //if(0 != modeId) // 目前只支持一种模型,信号对称拓延  'sym' or 'symh'  	Symmetric-padding (half-point): boundary value symmetric replication
        //	return -1;
        if( inLen < filterLen ) // inLen should lager than or equal extendLen, otherwise no extension
        	return -1;
    
        int i;
        int extendLen = filterLen - 1;
    
        if(0 == modeId) // 信号对称拓延
        {
            for(i=0; i<inLen; i++)
            {
            	out[extendLen+i] = in[i];
            }
            for(i=0; i<extendLen; i++)
            {
            	out[i]                     = out[2*extendLen - i - 1];       // 左边沿对称延拓
            	out[inLen + extendLen + i] = out[extendLen + inLen - i - 1]; // 右边沿对称延拓
            }
    
            return inLen + 2*extendLen;
        }
        else if(1 == modeId) // 信号周期拓延
        {
    		for( i = 0; i < extendLen; i++ )
    			out[i] = in[inLen-extendLen+i];
    		for ( i = 0; i < inLen; i++ )
    			out[extendLen+i] = in[i];
    
            return inLen + extendLen;
        }
    
    }

    1. /** 
    2.  * @brief 上采样  隔点插0 
    3.  * @param data 输入数据指针 
    4.  * @param n 输入数据长度 
    5.  * @param result 返回结果数组 
    6.  * @return 返回结果数组长度 
    7.  */  
    8. int Upsampling(double* data, int n, double result[])  
    9. {  
    10.     int i;  
    11.   
    12.     for( i = 0; i < n; i++ )  
    13.     {  
    14.         result[2*i] = data[i];  
    15.         result[2*i+1] = 0;  
    16.     }  
    17.   
    18.     return( 2*n );  
    19. }  
    /**
     * @brief 上采样  隔点插0
     * @param data 输入数据指针
     * @param n 输入数据长度
     * @param result 返回结果数组
     * @return 返回结果数组长度
     */
    int Upsampling(double* data, int n, double result[])
    {
    	int i;
    
    	for( i = 0; i < n; i++ )
    	{
    		result[2*i] = data[i];
    		result[2*i+1] = 0;
    	}
    
    	return( 2*n );
    }
    1. /** 
    2.  * @brief 下采样  隔点采样 
    3.  * @param data 输入数据指针 
    4.  * @param n 输入数据长度 
    5.  * @param result 返回结果数组 
    6.  * @return 返回结果数组长度 
    7.  */  
    8. int Downsampling(double* data, int n, double result[])  
    9. {  
    10.     int i, m;  
    11.   
    12.     m = n/2;  
    13.     for( i = 0; i < m; i++ )  
    14.         result[i] = data[i*2 + 1];  
    15.   
    16.     return( m );  
    17. }  
    /**
     * @brief 下采样  隔点采样
     * @param data 输入数据指针
     * @param n 输入数据长度
     * @param result 返回结果数组
     * @return 返回结果数组长度
     */
    int Downsampling(double* data, int n, double result[])
    {
    	int i, m;
    
    	m = n/2;
    	for( i = 0; i < m; i++ )
    		result[i] = data[i*2 + 1];
    
    	return( m );
    }


    1. /** 
    2.  * @brief 卷积运算 
    3.  * @param shapeId 卷积结果处理方式 
    4.  * @param double *inSignal, int signalLen, // 输入信号及其长度 
    5.  * @param double *inFilter, int filterLen, // 输入滤波器及其长度 
    6.  * @param double outConv[], int *convLen)   // 输出卷积结果及其长度 
    7.  * @return 
    8.  */  
    9. void Conv1(int shapeId,                  // 卷积结果处理方式  
    10.         double *inSignal, int signalLen, // 输入信号及其长度  
    11.         double *inFilter, int filterLen, // 输入滤波器及其长度  
    12.         double outConv[], int *convLen)   // 输出卷积结果及其长度  
    13. {  
    14.     if((NULL == inSignal)||(NULL == inFilter)||(NULL == outConv))  
    15.         return;  
    16.   
    17.     int n,k,kmin,kmax,p;  
    18.     if(0 == shapeId)      // 对于MATLAB conv(...,'shape')  -----full  
    19.     {  
    20.         *convLen = signalLen + filterLen - 1;  
    21.         for (n = 0; n < *convLen; n++)  
    22.         {  
    23.             outConv[n] = 0;  
    24.   
    25.             kmin = (n >= filterLen - 1) ? n - (filterLen - 1) : 0;  
    26.             kmax = (n < signalLen - 1) ? n : signalLen - 1;  
    27.   
    28.             for (k = kmin; k <= kmax; k++)  
    29.             {  
    30.                 outConv[n] += inSignal[k] * inFilter[n - k];  
    31.             }  
    32.         }  
    33.     }  
    34.     else if(1 == shapeId) // 对于MATLAB conv(...,'shape')  -----valid  
    35.     {  
    36.         *convLen = signalLen - filterLen + 1;  
    37.         for (n = filterLen - 1; n < signalLen; n++)  
    38.         {  
    39.             p = n - filterLen + 1;  
    40.             outConv[p] = 0;  
    41.   
    42.             kmin = (n >= filterLen - 1) ? n - (filterLen - 1) : 0;  
    43.             kmax = (n < signalLen - 1) ? n : signalLen - 1;  
    44.   
    45.             for (k = kmin; k <= kmax; k++)  
    46.             {  
    47.                 outConv[p] += inSignal[k] * inFilter[n - k];  
    48.             }  
    49.         }  
    50.     }  
    51.     else  
    52.         return ;  
    53.   
    54. }  
    /**
     * @brief 卷积运算
     * @param shapeId 卷积结果处理方式
     * @param double *inSignal, int signalLen, // 输入信号及其长度
     * @param double *inFilter, int filterLen, // 输入滤波器及其长度
     * @param double outConv[], int *convLen)   // 输出卷积结果及其长度
     * @return
     */
    void Conv1(int shapeId,                  // 卷积结果处理方式
    		double *inSignal, int signalLen, // 输入信号及其长度
    		double *inFilter, int filterLen, // 输入滤波器及其长度
    		double outConv[], int *convLen)   // 输出卷积结果及其长度
    {
        if((NULL == inSignal)||(NULL == inFilter)||(NULL == outConv))
            return;
    
        int n,k,kmin,kmax,p;
        if(0 == shapeId)      // 对于MATLAB conv(...,'shape')  -----full
        {
        	*convLen = signalLen + filterLen - 1;
        	for (n = 0; n < *convLen; n++)
        	{
        		outConv[n] = 0;
    
        	    kmin = (n >= filterLen - 1) ? n - (filterLen - 1) : 0;
        	    kmax = (n < signalLen - 1) ? n : signalLen - 1;
    
        	    for (k = kmin; k <= kmax; k++)
        	    {
        	    	outConv[n] += inSignal[k] * inFilter[n - k];
        	    }
        	}
        }
        else if(1 == shapeId) // 对于MATLAB conv(...,'shape')  -----valid
        {
        	*convLen = signalLen - filterLen + 1;
        	for (n = filterLen - 1; n < signalLen; n++)
        	{
        		p = n - filterLen + 1;
        		outConv[p] = 0;
    
        	    kmin = (n >= filterLen - 1) ? n - (filterLen - 1) : 0;
        	    kmax = (n < signalLen - 1) ? n : signalLen - 1;
    
        	    for (k = kmin; k <= kmax; k++)
        	    {
        	    	outConv[p] += inSignal[k] * inFilter[n - k];
        	    }
        	}
        }
        else
        	return ;
    
    }

    1. /** 
    2.  * @brief 小波变换之分解 
    3.  * @param sourceData 源数据 
    4.  * @param dataLen 源数据长度 
    5.  * @param db 过滤器类型 
    6.  * @param cA 分解后的近似部分序列-低频部分 
    7.  * @param cD 分解后的细节部分序列-高频部分 
    8.  * @return 正常则返回分解后序列的数据长度,错误则返回-1 
    9.  */  
    10. int Wavelet::Decomposition(double* sourceData, int dataLen, Filter db, double* cA, double* cD)  
    11. {  
    12.     if(dataLen < 2)  
    13.         return -1;  
    14.     if((NULL == sourceData)||(NULL == cA)||(NULL == cD))  
    15.         return -1;  
    16.   
    17.     m_db = db;  
    18.     int filterLen = m_db.length;  
    19.     int i, n;  
    20.     int decLen = (dataLen+filterLen-1)/2;  
    21.     int convLen = 0;  
    22.     double extendData[dataLen+2*filterLen-2];  
    23.     double convDataLow[dataLen+filterLen-1];  
    24.     double convDataHigh[dataLen+filterLen-1];  
    25.   
    26. /* 
    27. MATLAB上dwt函数的工作过程 
    28. 假设x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)],计算y=dwt(x,’db2’)。 
    29. 其计算过程主要由两个部分组成: 
    30. 1:边缘延拓,它主要由函数wextend完成。 
    31. 2:卷积运算,它主要由函数conv2完成。 
    32. 先看第一部分,仔细分析子程序部分,函数wextend的用法为y=wextend('1D','sym',x,3); 
    33. 这样得到的y=[ x(3) x(2) x(1) x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(7) x(6) x(5)] 
    34. 在看第二部分,仔细分析子程序部分,核心语句有z=conv2(y,Lo_D,'valid'); 
    35. 这里设Lo_D=[h(1) h(2) h(3) h(4)]。 
    36. 3:最后就是下采样,其下采样是按照式a = z(2:2:length(z))进行的,高频低频部分均如此,项数为floor((7+4-1)/2)。 
    37.  */  
    38.     // 1.边缘延拓  
    39.     SignalExtension(0, 0 , sourceData, dataLen, filterLen, extendData);  
    40.   
    41.     // 2.卷积运算  
    42.     Conv1(1, extendData, dataLen+2*filterLen-2, db.lowFilterDec, filterLen, convDataLow, &convLen);  
    43.     Conv1(1, extendData, dataLen+2*filterLen-2, db.highFilterDec, filterLen, convDataHigh, &convLen);  
    44.   
    45.     // 3.下采样  
    46.     Downsampling(convDataLow, dataLen + filterLen - 1, cA);  
    47.     Downsampling(convDataHigh, dataLen + filterLen - 1, cD);  
    48.   
    49.     return decLen;  
    50. }  
    /**
     * @brief 小波变换之分解
     * @param sourceData 源数据
     * @param dataLen 源数据长度
     * @param db 过滤器类型
     * @param cA 分解后的近似部分序列-低频部分
     * @param cD 分解后的细节部分序列-高频部分
     * @return 正常则返回分解后序列的数据长度,错误则返回-1
     */
    int Wavelet::Decomposition(double* sourceData, int dataLen, Filter db, double* cA, double* cD)
    {
        if(dataLen < 2)
            return -1;
        if((NULL == sourceData)||(NULL == cA)||(NULL == cD))
            return -1;
    
        m_db = db;
        int filterLen = m_db.length;
        int i, n;
        int decLen = (dataLen+filterLen-1)/2;
        int convLen = 0;
        double extendData[dataLen+2*filterLen-2];
        double convDataLow[dataLen+filterLen-1];
        double convDataHigh[dataLen+filterLen-1];
    
    /*
    MATLAB上dwt函数的工作过程
    假设x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)],计算y=dwt(x,’db2’)。
    其计算过程主要由两个部分组成:
    1:边缘延拓,它主要由函数wextend完成。
    2:卷积运算,它主要由函数conv2完成。
    先看第一部分,仔细分析子程序部分,函数wextend的用法为y=wextend('1D','sym',x,3);
    这样得到的y=[ x(3) x(2) x(1) x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(7) x(6) x(5)]
    在看第二部分,仔细分析子程序部分,核心语句有z=conv2(y,Lo_D,'valid');
    这里设Lo_D=[h(1) h(2) h(3) h(4)]。
    3:最后就是下采样,其下采样是按照式a = z(2:2:length(z))进行的,高频低频部分均如此,项数为floor((7+4-1)/2)。
     */
        // 1.边缘延拓
        SignalExtension(0, 0 , sourceData, dataLen, filterLen, extendData);
    
        // 2.卷积运算
        Conv1(1, extendData, dataLen+2*filterLen-2, db.lowFilterDec, filterLen, convDataLow, &convLen);
        Conv1(1, extendData, dataLen+2*filterLen-2, db.highFilterDec, filterLen, convDataHigh, &convLen);
    
        // 3.下采样
        Downsampling(convDataLow, dataLen + filterLen - 1, cA);
        Downsampling(convDataHigh, dataLen + filterLen - 1, cD);
    
        return decLen;
    }


    1. /** 
    2.  * @brief 小波变换之重构 
    3.  * @param cA 分解后的近似部分序列-低频部分 
    4.  * @param cD 分解后的细节部分序列-高频部分 
    5.  * @param cALength 输入数据长度 
    6.  * @param RecLength 输入重构后的原始数据长度 
    7.  * @param db 过滤器类型 
    8.  * @param recData 重构后输出的数据 
    9.  * @return 正常则返回重构数据长度,错误则返回-1 
    10.  */  
    11. int Wavelet::Reconstruction(double *cA, double *cD, int cALength, int RecLength, Filter db, double* recData)  
    12. {  
    13.     if((NULL == cA)||(NULL == cD)||(NULL == recData))  
    14.         return -1;  
    15.   
    16.     m_db = db;  
    17.     int filterLen = m_db.length;  
    18.   
    19.     int i,j;  
    20.     int n,k,p;  
    21.     int recLen = RecLength;  
    22.   
    23.     int convLen = 0;  
    24.     double convDataLow[recLen+filterLen-1];  
    25.     double convDataHigh[recLen+filterLen-1];  
    26.   
    27.     double cATemp[2*cALength];  
    28.     double cDTemp[2*cALength];  
    29.   
    30.     memset(convDataLow, 0, (recLen+filterLen-1)*sizeof(double)); // 清0  
    31.     memset(convDataHigh, 0, (recLen+filterLen-1)*sizeof(double)); // 清0  
    32.     memset(cATemp, 0, 2*cALength*sizeof(double)); // 清0  
    33.     memset(cDTemp, 0, 2*cALength*sizeof(double)); // 清0  
    34.   
    35.     // 1.隔点插0  
    36.     Upsampling(cA, cALength, cATemp);  
    37.     Upsampling(cD, cALength, cDTemp);  
    38.   
    39.     // 2.卷积运算  
    40.     Conv1(0, cATemp, 2*cALength-1, db.lowFilterRec, filterLen ,convDataLow, &convLen);  
    41.     convLen = 0;  
    42.     Conv1(0, cDTemp, 2*cALength-1, db.highFilterRec, filterLen ,convDataHigh, &convLen);  
    43.   
    44.     // 3.抽取结果及求和——实现类似MATLAB中的wkeep1(s,len,'c')的功能  
    45.     k = (convLen - recLen)/2;  
    46.     for(i=0; i<recLen; i++)  
    47.     {  
    48.         recData[i] = convDataLow[i + k] + convDataHigh[i + k];  
    49.     }  
    50.       
    51.     return recLen;  
    52. }  
    /**
     * @brief 小波变换之重构
     * @param cA 分解后的近似部分序列-低频部分
     * @param cD 分解后的细节部分序列-高频部分
     * @param cALength 输入数据长度
     * @param RecLength 输入重构后的原始数据长度
     * @param db 过滤器类型
     * @param recData 重构后输出的数据
     * @return 正常则返回重构数据长度,错误则返回-1
     */
    int Wavelet::Reconstruction(double *cA, double *cD, int cALength, int RecLength, Filter db, double* recData)
    {
        if((NULL == cA)||(NULL == cD)||(NULL == recData))
            return -1;
    
        m_db = db;
        int filterLen = m_db.length;
    
        int i,j;
        int n,k,p;
        int recLen = RecLength;
    
        int convLen = 0;
        double convDataLow[recLen+filterLen-1];
        double convDataHigh[recLen+filterLen-1];
    
        double cATemp[2*cALength];
        double cDTemp[2*cALength];
    
        memset(convDataLow, 0, (recLen+filterLen-1)*sizeof(double)); // 清0
        memset(convDataHigh, 0, (recLen+filterLen-1)*sizeof(double)); // 清0
        memset(cATemp, 0, 2*cALength*sizeof(double)); // 清0
        memset(cDTemp, 0, 2*cALength*sizeof(double)); // 清0
    
        // 1.隔点插0
        Upsampling(cA, cALength, cATemp);
        Upsampling(cD, cALength, cDTemp);
    
        // 2.卷积运算
        Conv1(0, cATemp, 2*cALength-1, db.lowFilterRec, filterLen ,convDataLow, &convLen);
        convLen = 0;
        Conv1(0, cDTemp, 2*cALength-1, db.highFilterRec, filterLen ,convDataHigh, &convLen);
    
        // 3.抽取结果及求和——实现类似MATLAB中的wkeep1(s,len,'c')的功能
        k = (convLen - recLen)/2;
        for(i=0; i<recLen; i++)
        {
        	recData[i] = convDataLow[i + k] + convDataHigh[i + k];
        }
    	
        return recLen;
    }

     

     

     

    //----转自  www.cnblogs.com/IDoIUnderstand/

     

    转载于:https://www.cnblogs.com/hjj-fighting/p/9767677.html

    展开全文
  • 离散序列的Mallat算法分解公式如下: 其中,H(n)、G(n)分别表示所选取的小波函数对应的低通和高通滤波器的抽头系数序列。 从Mallat算法的分解原理可知,分解后的序列就是原序列与滤波器序列的卷积再进行隔点...
  • 上次说到小波变换的知识体系,这篇博客就主要说小波变换里的连续信号的连续小波变换与离散小波变换。 连续信号的连续小波变换 话不多说,我们先放公式,如果你是第一次接触小波,你可会有点懵,但是不要怕,我希望...
  • 很容易拿傅里叶变换与小波变换对比着学习,但容易造成越比越混乱的现象,比如Matlab里fft函数所做的事就是离散傅里叶变换DFT,但Matlab里的dwt函数所做的事可不是离散小波变换的定义式,对于不想深入了解只想做应用...
  • 离散小波变换-DWT

    千次阅读 2014-10-13 10:47:29
    所谓的离散小波变换是指对尺度因子和时移因子
  • 离散小波变换DWT

    万次阅读 多人点赞 2017-04-28 23:38:31
    离散小波变换(Discrete Wavelet Transformation) 一、定义(摘自百度百科): 首先我们定义一些需要用到的信号及滤波器。x[n]:离散的输入信号,长度为N。g[n]:low pass filter低通滤波器,可以将输入信号的...
  • 这段代码实现的是MATLAB基于离散小波变换DWT)的语音和音频信号的数字水印代码,有GUI几乎每句都有详细的注释,附带一个录音的小程序,其中加入了两种干扰,一是低通滤波,二是白噪声干扰。另外还附赠解释小波变换...
  • 离散小波变换DWT)整理

    万次阅读 2016-04-20 18:05:02
    毕业设计涉及到对信号的去噪操作,而一维的离散小波变换可以运用在信号降噪中。因此对离散小波变换展开了调研,将内容整理如下。 离散小波变换(Discrete Wavelet Transformation) 一、定义(摘自百度百科): ...
  • 离散小波变换DWT)去噪

    千次阅读 2019-11-19 10:31:10
    小波去噪(DWT小波变换去噪优势方法 小波变换去噪优势 (1)能够很好的保留原有信号中所需的有用信号的峰值和突变部分,而用Fourior分析进行滤波时,由于有用信号集中在低频部分,而噪声集中在高频部分,通过低通...
  • 4层DWT离散小波变换

    2013-03-29 21:45:14
    使用97小波对图像进行4层小波变换,也可根据个人需要改变程序,对其进行任何层变换
  • 基于提升格式的离散小波变换比传统的基于卷积的运算量少,易于VLSI实现。本文提出了一种基于提升格式,高效实时实现JPEG2000中9/7双正交离散小波变换滤波器的VLSI结构设计方法。该方法所设计的结构,在保证同样的...
  • 【转】一维离散小波变换DWT)库,完全按matlab的wavelet toolbox 的API实现的 来源:http://hi.baidu.com/anatacollin/item/69fdab74ca7d045c0d0a07b4 一维离散小波变换DWT)库,完全按matlab的wavelet toolbox ...
  • 二维离散小波变换

    热门讨论 2012-10-31 22:38:49
    利用matlab程序实现二维离散小波变换,并对小波系数矩阵进行重构,进而在程序的编辑过程中理解二维离散小波变换和重构的原理和实现。 同时利用不同的小波和边缘延拓方法,对小波系数矩阵的能量、均值、方差、信噪比等...
  • 离散小波变换,用于压缩感知信号稀疏变换,是比较常用的稀疏基
  • 为什么会出先小波变换 窗口傅立叶变换(短时傅立叶变换)虽然可以部分定位时间,但由于窗口大小是固定的,只适用于频率波动小的平稳信号,不适用于频率波动大的非平稳信号。而小波变换可以根据频率的高低自动调节...
  • 离散小波变换

    2017-07-06 10:25:37
    离散小波变换的C程序 可运行
  • 【转】一维离散小波变换DWT)库,完全按matlab的wavelet toolbox 的API实现的 来源:http://hi.baidu.com/anatacollin/item/69fdab74ca7d045c0d0a07b4 一维离散小波变换DWT)库,完全按matlab的wavelet toolbox ...
  • 这是离散小波变换DWT)的MATLAB代码,利用db3离散小波变换对信号进行分解和重构,信号为一正弦叠加信号,采样点数2048,采样频率2000Hz,生成3个figure,分别是原始信号及其FFT,重构后的信号,还有重构信号的FFT

空空如也

1 2 3 4 5 ... 20
收藏数 1,916
精华内容 766
关键字:

离散小波变换