gsl库——小波变换

比眉伴天荒 2022-05-17 00:21 408阅读 0赞

### 1、定义 ###

本文将介绍gsl库的离散小波变换(Discrete Wavelet Transforms, DWTs)。本库包含针对一维或二维真实数据的小波处理。小波函数在头文件gsl\_wavelet.h和gsl\_wavelet2d.h中进行了声明。

连续小波变换及其逆变换公式如下：

![\\omega \\left ( s,\\tau \\right ) =\\int\_\{-\\infty \}^\{\\infty\} f(t) \* \\psi \_\{s,\\tau \}^\{\*\}(t)dt][omega _left _ s_tau _right _ _int_-_infty _infty_ f_t_ _ _psi _s_tau _t_dt]  
![f(t) =\\int\_\{0\}^\{\\infty\} ds \\int\_\{-\\infty \}^\{\\infty\} \\omega \\left ( s,\\tau \\right ) \* \\psi \_\{s,\\tau \} (t) d\\tau][f_t_ _int_0_infty_ ds _int_-_infty _infty_ _omega _left _ s_tau _right _ _ _psi _s_tau _ _t_ d_tau]

基函数![\\psi \_\{s,\\tau \}][psi _s_tau]通过一个单函数缩放(scaling)和平移(translation)得到，被称为母小波(mother wavelet)。  
小波变换的离散方法采用等距采样，采用固定的缩放和变换步长![(s,\\tau )][s_tau]。频率和时间轴采用二分采样的方法，缩放倍率为![2^\{j\}][2_j]倍，参数为j。  
函数族![\\\{\\psi \_\{s,n \}\\\}][psi _s_n]的结果族组成了平方可积(square-integrable)信号的标准正交基(orthonormal basis)。离散小波变换是一种![O\\left ( N \\right )][O_left _ N _right]复杂度算法，也被称为快速小波变换(fast wavelet transform)。

### 2、初始化 ###

**gsl\_wavelet**

这个结构包含了定义小波以及任何相关偏移参数(offset parameter)的滤波系数(filter coefficient)。

**gsl\_wavelet \* gsl\_wavelet\_alloc(const gsl\_wavelet\_type \* T, size\_t k)**  
此函数分配并初始化了一个小波对象T，参数k选择小波族里的特定类型，内存不足或不支持类型返回空指针。

**gsl\_wavelet\_type  
gsl\_wavelet\_daubechies  
gsl\_wavelet\_daubechies\_centered**  
这是消失矩(vanishing moments)k＝2的最大相位的Daubechies小波族。实现小波是k＝4，6，...，20，k为偶数

**gsl\_wavelet\_haar  
gsl\_wavelet\_haar\_centered**  
这是Haar小波。Haar小波唯一有效选择是k＝2。

**gsl\_wavelet\_bspline  
gsl\_wavelet\_bspline\_centered**  
阶数为![(i,j)][i_j]的双正交B样条(biorthogonal B-spline) 小波族(wavelet family) 。 实现值![k = 100 \* i+j][k _ 100 _ i_j]包括103, 105, 202, 204, 206, 208, 301, 303, 305 307, 309.  
小波的中心形式将边缘上的各种子带的系数对齐。从而导致在相平面中小波变换系数的可视化更容易理解。

**const char \* gsl\_wavelet\_name(const gsl\_wavelet \* w)**  
此函数返回了一个指向小波族成员的指针w。

**void gsl\_wavelet\_free(gsl\_wavelet \* w)**  
此函数释放了小波对象w。

**gsl\_wavelet\_workspace**  
此结构包含与输入数据相同大小的暂存空间(scratch space)，并用于在转换过程中保存中间结果。

**gsl\_wavelet\_workspace \* gsl\_wavelet\_workspace\_alloc(size\_t n)**  
该函数分配离散小波变换的工作空间。对n个元素执行一维变换，必须提供一个大小为n的工作空间。关于n×n矩阵的二维变换分配大小为n的工作空间已足够，因为变换在单独的行和列上运行。内存不足返回空指针。

**void gsl\_wavelet\_workspace\_free(gsl\_wavelet\_workspace \* work)**  
该函数释放为work分配的工作空间。

### 3、变换函数 ###

这部分描述了执行离散小波变换的实际函数。注意，变换使用周期边界条件(periodic boundary conditions)。如果信号在样本长度上不是周期性的，那么在变换的每个级别的开始和结束时都会出现杂散系数(spurious coefficients)。

**3.1一维小波变换**

**int gsl\_wavelet\_transform(const gsl\_wavelet \* w, double \* data, size\_t stride, size\_t n, gsl\_wavelet\_direction dir, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet\_transform\_forward(const gsl\_wavelet \* w, double \* data, size\_t stride, size\_t n, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet\_transform\_inverse(const gsl\_wavelet \* w, double \* data, size\_t stride, size\_t n, gsl\_wavelet\_workspace \* work)**

这些函数计算了长度为*n、*步长为*stride*的正逆离散小波变换。变换长度n限于2的幂次。对于函数的变换版本，参数dir可以是forward(+1)或backward(-1)。必须提供长度为n的工作空间。  
对正变换，原始数组的元素被离散小波变换以![f\_\{i\}\\rightarrow \\omega \_\{j,k\}][f_i_rightarrow _omega _j_k]填充三角形存储布局(packed triangular storage layout)取代，j是迭代层次![j=0...J-1][j_0...J-1]的指数，k是每层系数![k=0...2^\{j\}-1][k_0...2_j_-1]的指数。层次总数是![J=log\_\{2\}(n)][J_log_2_n]。输出数据有一下形式：

![\\left ( s\_\{-1,0\},d\_\{0,0\},d\_\{1,0\}, ...,d\_\{j,k\},...,d\_\{J-1,2^\{J-1\}-1\}\\right )][left _ s_-1_0_d_0_0_d_1_0_ ..._d_j_k_..._d_J-1_2_J-1_-1_right]  
第一个参数![s\_\{-1,0\}][s_-1_0]是光滑因子，每层j的细节因子![d\_\{j,k\}][d_j_k]。逆变换逆换这些系数来获取原始数据。  
这些函数成功运行完成则返回状态参数**GSL\_SUCCESS**。如果n不是2的幂次或内存不足则返回**GSL\_EINVAL**。

**3.2二维小波变换**

**int gsl\_wavelet2d\_transform(const gsl\_wavelet \* w, double \* data, size\_t tda, size\_t size1, size\_t size2, gsl\_wavelet\_direction dir, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet2d\_transform\_forward(const gsl\_wavelet \* w, double \* data, size\_t tda, size\_t size1, size\_t size2, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet2d\_transform\_inverse(const gsl\_wavelet \* w, double \* data, size\_t tda, size\_t size1, size\_t size2, gsl\_wavelet\_workspace \* work)**

**int gsl\_wavelet2d\_transform\_matrix(const gsl\_wavelet \* w, gsl\_matrix \* m, gsl\_wavelet\_direction dir, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet2d\_transform\_matrix\_forward(const gsl\_wavelet \* w, gsl\_matrix \* m, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet2d\_transform\_matrix\_inverse(const gsl\_wavelet \* w, gsl\_matrix \* m, gsl\_wavelet\_workspace \* work)**

**int gsl\_wavelet2d\_nstransform(const gsl\_wavelet \* w, double \* data, size\_t tda, size\_t size1, size\_t size2, gsl\_wavelet\_direction dir, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet2d\_nstransform\_forward(const gsl\_wavelet \* w, double \* data, size\_t tda, size\_t size1, size\_t size2, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet2d\_nstransform\_inverse(const gsl\_wavelet \* w, double \* data, size\_t tda, size\_t size1, size\_t size2, gsl\_wavelet\_workspace \* work)**

**int gsl\_wavelet2d\_nstransform\_matrix(const gsl\_wavelet \* w, gsl\_matrix \* m, gsl\_wavelet\_direction dir, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet2d\_nstransform\_matrix\_forward(const gsl\_wavelet \* w, gsl\_matrix \* m, gsl\_wavelet\_workspace \* work)  
int gsl\_wavelet2d\_nstransform\_matrix\_inverse(const gsl\_wavelet \* w, gsl\_matrix \* m, gsl\_wavelet\_workspace \* work)**

### 4、代码示例 ###

以下程序使用了一维小波变换函数。本例对256长度的输入信号近似处理，通过小波变换，保留20个最大的成分，将剩余成分置零。

#include <stdio.h>
    #include <math.h>
    #include <gsl/gsl_sort.h>
    #include <gsl/gsl_wavelet.h>
    int
    main (int argc, char **argv)
    {
        (void)(argc); /* avoid unused parameter warning */
        int i, n = 256, nc = 20;
        double *orig_data = malloc (n * sizeof (double));
        double *data = malloc (n * sizeof (double));
        double *abscoeff = malloc (n * sizeof (double));
        size_t *p = malloc (n * sizeof (size_t));
        FILE * f;
        gsl_wavelet *w;
        gsl_wavelet_workspace *work;
        w = gsl_wavelet_alloc (gsl_wavelet_daubechies, 4);
        work = gsl_wavelet_workspace_alloc (n);
        f = fopen (argv[1], "r");
        for (i = 0; i < n; i++)
        {
            fscanf (f, "%lg", &orig_data[i]);
            data[i] = orig_data[i];
        }
        fclose (f);
        gsl_wavelet_transform_forward (w, data, 1, n, work);
        for (i = 0; i < n; i++)
        {
            abscoeff[i] = fabs (data[i]);
        }
        gsl_sort_index (p, abscoeff, 1, n);
        for (i = 0; (i + nc) < n; i++)
            data[p[i]] = 0;
        gsl_wavelet_transform_inverse (w, data, 1, n, work);
        for (i = 0; i < n; i++)
        {
            printf ("%g %g\n", orig_data[i], data[i]);
        }
        gsl_wavelet_free (w);
        gsl_wavelet_workspace_free (work);
        free (data);
        free (orig_data);
        free (abscoeff);
        free (p);
        return 0;
    }

### 5、小波降噪 ###

void waveletDenoise(double *pdData, int iDataLen)
    {
        gsl_wavelet *w = NULL;
        gsl_wavelet_workspace *work = NULL;
        double *abscoeff = (double *)malloc(iDataLen * sizeof (double));
        int i=0, iFrom=0, iTo=0, j=0;
        int iLevel = binary_logn(iDataLen);
        double dSign = 0.0, dThres = 0.0, dCoefDif = 0.0, dCoeffstd = 0.0;
    
    
        if (NULL == abscoeff)
        {
            return ;
        }
    
        w    = gsl_wavelet_alloc(gsl_wavelet_daubechies, 10);
        work = gsl_wavelet_workspace_alloc(iDataLen);
        if ((NULL != w) &&
            (NULL != work))
        {
            gsl_wavelet_transform_forward(w, pdData, 1, iDataLen, work);
            for (i = 0; i < iDataLen; i++)
            {
                abscoeff[i] = fabs(pdData[i]);
            }
    
            for (i=8; i<iLevel; i++)
            {
                iFrom  = pow(2, i);
                iTo    = pow(2, i+1) - 1;
    
                dCoeffstd = calcMthMax(abscoeff, iFrom, iTo, (iTo - iFrom + 1) / 2) / 0.6745;
    
                dThres = 0.3936 + 0.1829*(i);
                dThres *= dCoeffstd;
                
                for (j=iFrom; j<iTo; j++)
                {
                    if (pdData[j] < 0)
                    {
                        dSign = -1.0;
                    }
                    else if (pdData[j] > 0)
                    {
                        dSign = 1.0;
                    }
    
                    dCoefDif  = abscoeff[j] - dThres;
                    dCoefDif  = (dCoefDif + fabs(dCoefDif)) / 2;
                    pdData[j] = dSign * dCoefDif;
                }
            }
    
            gsl_wavelet_transform_inverse (w, pdData, 1, iDataLen, work);                
        }
    
        gsl_wavelet_free (w);
        gsl_wavelet_workspace_free (work);
       
    
        if (NULL != abscoeff)
        {
            free(abscoeff);
            abscoeff = NULL;
        }
        
        return ;
    }

[omega _left _ s_tau _right _ _int_-_infty _infty_ f_t_ _ _psi _s_tau _t_dt]: https://private.codecogs.com/gif.latex?%5Comega%20%5Cleft%20%28%20s%2C%5Ctau%20%5Cright%20%29%20%3D%5Cint_%7B-%5Cinfty%20%7D%5E%7B%5Cinfty%7D%20f%28t%29%20*%20%5Cpsi%20_%7Bs%2C%5Ctau%20%7D%5E%7B*%7D%28t%29dt
[f_t_ _int_0_infty_ ds _int_-_infty _infty_ _omega _left _ s_tau _right _ _ _psi _s_tau _ _t_ d_tau]: https://private.codecogs.com/gif.latex?f%28t%29%20%3D%5Cint_%7B0%7D%5E%7B%5Cinfty%7D%20ds%20%5Cint_%7B-%5Cinfty%20%7D%5E%7B%5Cinfty%7D%20%5Comega%20%5Cleft%20%28%20s%2C%5Ctau%20%5Cright%20%29%20*%20%5Cpsi%20_%7Bs%2C%5Ctau%20%7D%20%28t%29%20d%5Ctau
[psi _s_tau]: https://private.codecogs.com/gif.latex?%5Cpsi%20_%7Bs%2C%5Ctau%20%7D
[s_tau]: https://private.codecogs.com/gif.latex?%28s%2C%5Ctau%20%29
[2_j]: https://private.codecogs.com/gif.latex?2%5E%7Bj%7D
[psi _s_n]: https://private.codecogs.com/gif.latex?%5C%7B%5Cpsi%20_%7Bs%2Cn%20%7D%5C%7D
[O_left _ N _right]: https://private.codecogs.com/gif.latex?O%5Cleft%20%28%20N%20%5Cright%20%29
[i_j]: https://private.codecogs.com/gif.latex?%28i%2Cj%29
[k _ 100 _ i_j]: https://private.codecogs.com/gif.latex?k%20%3D%20100%20*%20i&plus;j
[f_i_rightarrow _omega _j_k]: https://private.codecogs.com/gif.latex?f_%7Bi%7D%5Crightarrow%20%5Comega%20_%7Bj%2Ck%7D
[j_0...J-1]: https://private.codecogs.com/gif.latex?j%3D0...J-1
[k_0...2_j_-1]: https://private.codecogs.com/gif.latex?k%3D0...2%5E%7Bj%7D-1
[J_log_2_n]: https://private.codecogs.com/gif.latex?J%3Dlog_%7B2%7D%28n%29
[left _ s_-1_0_d_0_0_d_1_0_ ..._d_j_k_..._d_J-1_2_J-1_-1_right]: https://private.codecogs.com/gif.latex?%5Cleft%20%28%20s_%7B-1%2C0%7D%2Cd_%7B0%2C0%7D%2Cd_%7B1%2C0%7D%2C%20...%2Cd_%7Bj%2Ck%7D%2C...%2Cd_%7BJ-1%2C2%5E%7BJ-1%7D-1%7D%5Cright%20%29
[s_-1_0]: https://private.codecogs.com/gif.latex?s_%7B-1%2C0%7D
[d_j_k]: https://private.codecogs.com/gif.latex?d_%7Bj%2Ck%7D