C语言实现一维信号小波阙值去噪

实现小波阙值去噪分三个步骤：

1.将源信号分解成系数；

2.根据分解的系数获取阙值；

3.根据阙值对系数进行过滤；

4.将过滤的系数还原成信号，最终的结果就是将源信号过滤后的结果。

我实现的代码主要是参阅了这篇文章：https://blog.csdn.net/EbowTang/article/details/40481393

想了解小波阙值去噪的概念或者更多信息都可以去这篇文章点赞参观。

接下来我就附上源码：

编译命令例子：gcc src.c -o src.o -lm

#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
#include <math.h>

#define LENGTH 10240//信号长度
#define READFILE "1_laser.txt" /*源信号文件*/

typedef struct tagCL1D
{//本结构体用于接收一维图像分解信息，避免参数的大量传送
       int Scale;//分解尺度
       int dbn;//db小波的编号
       int allSize;//总数据长度
       int * msgLen;
}msgCL1D;

msgCL1D m_msgCL1D;

//小波滤波器
struct tagWaveFilter
{
       int filterLen;
       double *Lo_D;
       double *Hi_D;
       double *Lo_R;
       double *Hi_R;
};

double db6_Lo_D[12] = { -0.00107730108499558   ,0.00477725751101065,   0.000553842200993802, - 0.0315820393180312,   0.0275228655300163,0.0975016055870794 ,- 0.129766867567096 ,- 0.226264693965169,   0.315250351709243,   0.751133908021578,   0.494623890398385,   0.111540743350080};
double db6_Hi_D[12] = { -0.111540743350080,   0.494623890398385, - 0.751133908021578,   0.315250351709243,   0.226264693965169,- 0.129766867567096, - 0.0975016055870794,   0.0275228655300163,   0.0315820393180312,   0.000553842200993802, - 0.00477725751101065 ,- 0.00107730108499558};
double db6_Lo_R[12] = { 0.111540743350080,   0.494623890398385,   0.751133908021578,   0.315250351709243, - 0.226264693965169 ,- 0.129766867567096,0.0975016055870794   ,0.0275228655300163, - 0.0315820393180312,   0.000553842200993802,   0.00477725751101065 ,- 0.00107730108499558};
double db6_Hi_R[12] = { -0.00107730108499558 ,- 0.00477725751101065,   0.000553842200993802,   0.0315820393180312,   0.0275228655300163, - 0.0975016055870794,- 0.129766867567096,0.226264693965169,   0.315250351709243, - 0.751133908021578,   0.494623890398385, - 0.111540743350080};
double db3_Lo_D[6] = { 0.0352262918821007, -0.0854412738822415, -0.135011020010391, 0.459877502119331, 0.806891509313339, 0.332670552950957 };
double db3_Hi_D[6] = { -0.332670552950957, 0.806891509313339, -0.459877502119331, -0.135011020010391, 0.0854412738822415, 0.0352262918821007 };
double db3_Lo_R[6] = { 0.332670552950957, 0.806891509313339, 0.459877502119331, -0.135011020010391, -0.0854412738822415, 0.0352262918821007 };
double db3_Hi_R[6] = { 0.0352262918821007, 0.0854412738822415, -0.135011020010391, -0.459877502119331, 0.806891509313339, -0.332670552950957 };
double db10_Lo_D[20] = { -0.0000, 0.0001, -0.0001, -0.0007, 0.0020, 0.0014, -0.0107,0.0036, 0.0332, -0.0295, -0.0714, 0.0931, 0.1274, -0.1959, -0.2498, 0.2812,0.6885, 0.5272, 0.1882, 0.0267 };
double db10_Hi_D[20] = { -0.0267, 0.1882, -0.5272, 0.6885, -0.2812, -0.2498, 0.1959,0.1274, -0.0931, -0.0714, 0.0295, 0.0332, -0.0036, -0.0107, -0.0014, 0.0020,0.0007, -0.0001, -0.0001, -0.0000 };
double db10_Lo_R[20] = { 0.0267, 0.1882, 0.5272, 0.6885, 0.2812, -0.2498, -0.1959,0.1274, 0.0931, -0.0714, -0.0295, 0.0332, 0.0036, -0.0107, 0.0014, 0.0020,-0.0007, -0.0001, 0.0001, -0.0000 };
double db10_Hi_R[20] = { -0.0000, -0.0001, -0.0001, 0.0007, 0.0020, -0.0014, -0.0107,-0.0036, 0.0332, 0.0295, -0.0714, -0.0931, 0.1274, 0.1959, -0.2498, -0.2812,0.6885, -0.5272, 0.1882, -0.0267 };
double db1_Lo_D[2] = { 0.7071, 0.7071 };
double db1_Hi_D[2] = { -0.7071, 0.7071 };
double db1_Lo_R[2] = { 0.7071, 0.7071 };
double db1_Hi_R[2] = { 0.7071, -0.7071 };
struct tagWaveFilter m_dbFilter;

// 初始化一维小波分解信息，从中保存未来需要的信息
int InitDecInfo(const int signalLen,const int decScale,const int decdbn)
/*const int signalLen,源信号长度;const int decScale,分解尺度;const int decdbn,db滤波器的编号*/
{

if (signalLen < m_dbFilter.filterLen - 1)
   {
       printf("错误信息：滤波器长度大于信号!\n");
       return -1;
   }

int srcLen = signalLen;
   m_msgCL1D.dbn = decdbn;
   m_msgCL1D.Scale = decScale;
   //m_msgCL1D.msgLen.resize(decScale + 2);
   m_msgCL1D.msgLen=(int *)calloc(decScale + 2,sizeof(int));
   if(m_msgCL1D.msgLen==NULL)
   {
       printf("%s[%d]:calloc error\n",__FUNCTION__,__LINE__);
       return -1;
   }
   m_msgCL1D.msgLen[0] = srcLen;
   for (int i = 1; i <= decScale; i++)
   {
       int exLen = (srcLen + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度
       srcLen = exLen;
       m_msgCL1D.msgLen[i] = srcLen;
   }
   m_msgCL1D.msgLen[decScale + 1] = srcLen;

for (int i = 1; i < decScale + 2; i++)
       m_msgCL1D.allSize += m_msgCL1D.msgLen[i];

return 0;
}

// 一维信号的小波分解
/*   double *pSrcData,分解的源信号;int srcLen,源信号的长度;double *pDstCeof分解出来的,本函数将返回此长度*/
int DWT(double *pSrcData,int srcLen,double *pDstCeof)
{
   //本程序禁止出现这种情况，否则数据出错（对称拓延长度为filterLen-1，如果大于了signalLen将越界）
   if (srcLen < m_dbFilter.filterLen - 1)
   {   //实际上信号的长度可以是任意的（matlab顺序：信号拓延-》卷积-》下采样），
       //但是本程序为了算法速度，写法上不允许
       printf("错误信息：滤波器长度大于信号!\n");
       exit(1);
   }
   int exLen = (srcLen + m_dbFilter.filterLen - 1) / 2;//对称拓延后系数的长度
   int k = 0;
   double tmp = 0.0;
   for (int i = 0; i < exLen; i++)
   {

pDstCeof[i] = 0.0;
       pDstCeof[i + exLen] = 0.0;
       for (int j = 0; j < m_dbFilter.filterLen; j++)
       {
           k = 2 * i - j + 1;
           //信号边沿对称延拓
           if ((k<0) && (k >= -m_dbFilter.filterLen + 1))//左边沿拓延
               tmp = pSrcData[-k - 1];
           else if ((k >= 0) && (k <= srcLen - 1))//保持不变
               tmp = pSrcData[k];
           else if ((k>srcLen - 1) && (k <= (srcLen + m_dbFilter.filterLen - 2)))//右边沿拓延
               tmp = pSrcData[2 * srcLen - k - 1];
           else
               tmp = 0.0;
           pDstCeof[i] += m_dbFilter.Lo_D[j] * tmp;
           pDstCeof[i + exLen] += m_dbFilter.Hi_D[j] * tmp;
       }
   }
   return 2 * exLen;
}

// 计算阈值
/*   double *pDetCoef,细节系数（应该是第一级的细节系数）;int detLen,此段细节系数的长度*/
double getThr(double *pDetCoef,int detLen)
{
   double thr = 0.0;
   double sigma = 0.0;
   int i,j;
   double tmp;
   for ( i = 0; i < detLen; i++)
       pDetCoef[i] = fabs(pDetCoef[i]);

//std::sort(pDetCoef, pDetCoef + detLen);
   for ( i = 0; i < detLen-1; i++)
   {
       for(j=0;j<detLen-i-1;j++)
       {
           if(pDetCoef[j]>pDetCoef[j+1])
           {
               tmp=pDetCoef[j+1];
               pDetCoef[j+1]=pDetCoef[j];
               pDetCoef[j]=tmp;
           }
       }
   }

if (detLen % 2 == 0 && detLen >= 2)
       sigma = (pDetCoef[detLen / 2-1] + pDetCoef[detLen / 2]) / 2 / 0.6745;
   else
       sigma = pDetCoef[detLen / 2] / 0.6745;

double N = m_msgCL1D.msgLen[0];
   thr = sigma *sqrt(2.0*log(N));

return thr;
}

// 将系数阈值处理，一维二维均适用
/*   double *pDstCoef,细节系数（应该是除近似系数外的所有的细节系数）;double thr,阈值;const int allsize,分解出来的系数的总长度（非）;const int gap,跳过最后一层的近似系数*/
void Wthresh(double *pDstCoef,double thr,const int allsize,const int gap,int ish)
{
   if (ish)//硬阈值
   {
       for (int i = gap; i < allsize; i++)
       {
           if (fabs(pDstCoef[i]) < thr)
               pDstCoef[i] = 0.0;
       }
   }
   else//软阈值
   {
       for (int i = gap; i < allsize; i++)
       {
           if (fabs(pDstCoef[i]) < thr)
           {
               pDstCoef[i] = 0.0;
           }
           else
           {
               if (pDstCoef[i] < 0.0)
                   pDstCoef[i] = thr - fabs(pDstCoef[i]);
               else
                   pDstCoef[i] = fabs(pDstCoef[i]) - thr;

           }
       }
   }
}
// 将系数阈值处理，一维二维均适用
/*   double *pDstCoef,细节系数（应该是除近似系数外的所有的细节系数）;double thr,阈值;const int allsize,分解出来的系数的总长度（非）;const int gap,跳过最后一层的近似系数*/
void myWthresh(double *pDstCoef,double thr,const int num,int ish)
{
   if (ish)//硬阈值
   {
       for (int i = 0; i < num; i++)
       {
           if (fabs(pDstCoef[i]) < thr)
               pDstCoef[i] = 0.0;
       }
   }
   else//软阈值
   {
       for (int i = 0; i < num; i++)
       {
           if (fabs(pDstCoef[i]) < thr)
           {
               pDstCoef[i] = 0.0;
           }
           else
           {
               if (pDstCoef[i] < 0.0)
                   pDstCoef[i] = thr - fabs(pDstCoef[i]);
               else
                   pDstCoef[i] = fabs(pDstCoef[i]) - thr;

           }
       }
   }
}
// 一维小波反变换，重构出源信号
/*   double *pSrcCoef,源分解系数;int dstLen,重构出来的系数的长度;double *pDstData重构出来的系数*/
void IDWT(double *pSrcCoef,int dstLen,double *pDstData)
{
   int p = 0;
   int caLen = (dstLen + m_dbFilter.filterLen - 1) / 2;
   for (int i = 0; i < dstLen; i++)
   {
       pDstData[i] = 0.0;
       for (int j = 0; j < caLen; j++)
       {
           p = i - 2 * j + m_dbFilter.filterLen - 2;
           //信号重构
           if ((p >= 0) && (p<m_dbFilter.filterLen))
               pDstData[i] += m_dbFilter.Lo_R[p] * pSrcCoef[j] + m_dbFilter.Hi_R[p] * pSrcCoef[j + caLen];
       }
   }
}

上面这一段是头文件，以及main函数需要用到的代码。

下面是我写的一些main函数用于不同的测试（例子就写一个）：

1.将小波阙值去噪封装成函数generatewave，直接调该函数即可，代码都可以拿去直接用。

/*dstsrc 为源信息，n为源信息的个数，dstsrc 为最后的处理结果,在调用函数中dstsrc要开辟n个double类型的空间并且要初始化，若为空则会报错*/
int generatewave(double *srcmes,int n,double *dstsrc)
{
   double *dst=NULL;
   int dstnum;
   double *dsttmp=NULL;
   double thr;

   if(srcmes==NULL)
   {
       printf("%s[%d]:srcmes is null\n",__FUNCTION__,__LINE__);
       return -1;
   }
   m_dbFilter.filterLen=10;
   m_dbFilter.Lo_D=db10_Lo_D;
   m_dbFilter.Hi_D=db10_Hi_D;
   m_dbFilter.Lo_R=db10_Lo_R;
   m_dbFilter.Hi_R=db10_Hi_R;

   InitDecInfo(n,10,10);

   dst = (double *)calloc(m_msgCL1D.allSize,sizeof(double));
   dstnum=DWT(srcmes, n,dst);

   dsttmp = (double *)calloc(dstnum,sizeof(double));
   for (int i = 0; i < dstnum; i++)
       dsttmp[i]=dst[i];

   thr=getThr(dsttmp,dstnum);

   /*for (int i = 0; i < 10; i++)
       printf("%lf ",dst[i]);
   printf("\n");*/

   myWthresh(dst,thr,dstnum,0);

   //dstsrc = (double *)calloc(n,sizeof(double));
   if(dstsrc==NULL)
   {
       printf("%s[%d]:dstsrc is null\n",__FUNCTION__,__LINE__);
       return -1;
   }
   IDWT(dst, n, dstsrc);
   for (int i = 0; i < n; i++)
       printf("%lf ",dstsrc[i]);
   printf("\n");
   if(dst)
       free(dst);
   if(dsttmp)
       free(dsttmp);
   return 0;
}

int main()
{
       //读取输入信号
   char s[64];
FILE *fp;
double *data=NULL;
int n=0;
       double *dstsrc=NULL;
fp=fopen(READFILE,"r");
if(fp==NULL) //如果读取失败
{
printf("错误！找不到要读取的文件\"%s\"\n",READFILE);
exit(1);//中止程序
}
bzero(s,64);
    while( fgets(s, 32, fp) != NULL )//读取长度n要设置得长一点，要保证读到回车符，这样指针才会定位到下一行？回车符返回的是零值？是，非数字字符经过atoi变换都应该返回零值
{
       data=(double *)realloc(data,(n+1)*sizeof(double));
data[n] = atof(s);
n++;
printf("[%d]:%s,%lf\n",n,s,data[n-1]);
bzero(s,64);
}
fclose(fp);
   printf("ok====\n");

   dstsrc = (double *)calloc(n,sizeof(double));
   if(0!=generatewave( data, n, dstsrc))
   {
       return -1;
   }

   fp=fopen("data_output.txt","w");
   if(fp==NULL)
   {
       if(dstsrc)
           free(dstsrc);
       if(data)
           free(data);
       printf("%s[%d]:\n",__FUNCTION__,__LINE__);
       return -1;
   }
   if(dstsrc==NULL)
   {
       printf("%s[%d]:dstsrc is null\n",__FUNCTION__,__LINE__);
   }
//打印一维小波变换后的结果
for(int i = 0; i < n; i++)
{
   printf("%s[%d]:\n",__FUNCTION__,__LINE__);
fprintf(fp,"%lf\n", dstsrc[i]);
}
   fclose(fp);
   for (int i = 0; i < n; i++)
           printf("%lf ",dstsrc[i]);
   printf("\n");
   printf("over===\n");

   if(dstsrc)
       free(dstsrc);
   if(data)
       free(data);

return 0;
}

这里附上我处理的源码信息共9997行：

C语言实现一维信号小波阙值去噪相关推荐

一维信号小波阈值去噪 c语言,一维信号小波阈值去噪
描述 1.小波阈值处理基本理论所谓阈值去噪简而言之就是对信号进行分解,然后对分解后的系数进行阈值处理,最后重构得到去噪信号.该算法其主要理论依据是:小波变换具有很强的去数据相关性,它能够使信号的能量 ...
一维信号小波阈值去噪
转载:http://blog.csdn.net/ebowtang/article/details/40481393 一,小波阈值去噪基本理论本博文根据小波的分解与重构原理,实现了基于硬阈 ...
【心电信号】基于matlab小波阙值心电信号去噪【含Matlab源码 2188期】
⛄一.小波阈值法去噪概述电能质量扰动信号的噪声大多以高斯白噪声的形式存在,利用小波变换对信号进行多分辨率分解,由于小波变换具有去除数据相关性的特点,故可以将有用信号与噪声的能量分离开来.信号中有效的 ...
【数字信号去噪】基于matlab小波阙值数字信号去噪和求信噪比【含Matlab源码 2191期】
⛄一.小波阈值法去噪概述电能质量扰动信号的噪声大多以高斯白噪声的形式存在,利用小波变换对信号进行多分辨率分解,由于小波变换具有去除数据相关性的特点,故可以将有用信号与噪声的能量分离开来.信号中有效的 ...
【图像去噪】基于自适应小波阙值算法实现图像去噪附matlab代码
1 内容介绍基于 Donoho经典小波阈值去除图像噪声基本思路,分析常用硬阈值法和软阈值法在图像去噪中的缺陷.针对这些缺陷,提出一种改进的阈值去噪法,该方法不仅可克服硬阈值不连续的缺点,还能够有效解 ...
一维信号小波去噪算法C语言,[转载]一维小波分解与去噪重构
对随机一维信号实现多尺度小波分解,长度M=256,层数N=3. 解:一维随机信号是用nelec函数作为信号源,产生0~256范围内的信号. 实现小波按层分解的函数是: [C,L]=wavedec(s, ...
一维信号小波去噪原理及python实现示例
信号去噪是经常用到的信号预处理过程,以达到在保留原有信号真是信息的基础上尽可能低降低或者消除噪声,获得更高质量的信号,从而为下一步的处理奠定基础. 去噪方法可分为时域方法与频域方法.时域方法是指直接在 ...
matlab非平稳信号小波和FFT去噪
文章目录目录文章目录一.非平稳信号二.FFT 和小波去噪实例完整代码一.非平稳信号在实际的工程应用中,大多数信号可能包含着许多尖峰或突变,而且噪声信号也并不是平衡的白噪声.对这种信号进行 ...
舰船辐射噪声 matlab,基于MATLAB的舰船辐射噪声信号小波消噪处理
基于MATLAB的舰船辐射噪声信号小波消噪处理针对舰船辐射噪声信号的特点提出了小波消噪的方法,对小波消噪理论作了简要的阐述,并设计了一种消噪方案,最后利用MATLAB,在计算机 (本文共4页) 阅读 ...

C语言实现一维信号小波阙值去噪

C语言实现一维信号小波阙值去噪相关推荐

最新文章

热门文章