gsl库——小波变换

最新推荐文章于 2025-03-19 17:17:02 发布

danxibaoxxx

最新推荐文章于 2025-03-19 17:17:02 发布

阅读量3.4k

点赞数

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$
$f(t) =\int_{0}^{\infty} ds \int_{-\infty }^{\infty} \omega \left ( s,\tau \right ) * \psi _{s,\tau } (t) d\tau$

基函数 $\psi _{s,\tau }$ 通过一个单函数缩放(scaling)和平移(translation)得到，被称为母小波(mother wavelet)。
小波变换的离散方法采用等距采样，采用固定的缩放和变换步长 $(s,\tau )$ 。频率和时间轴采用二分采样的方法，缩放倍率为 $2^{j}$ 倍，参数为j。
函数族 $\{\psi _{s,n }\}$ 的结果族组成了平方可积(square-integrable)信号的标准正交基(orthonormal basis)。离散小波变换是一种 $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)$ 的双正交B样条(biorthogonal B-spline) 小波族(wavelet family) 。实现值 $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}$ 填充三角形存储布局(packed triangular storage layout)取代，j是迭代层次 $j=0...J-1$ 的指数，k是每层系数 $k=0...2^{j}-1$ 的指数。层次总数是 $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 )$
第一个参数 $s_{-1,0}$ 是光滑因子，每层j的细节因子 $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 ;
}