首页 | 新闻 | 新品 | 文库 | 方案 | 视频 | 下载 | 商城 | 开发板 | 数据中心 | 座谈新版 | 培训 | 工具 | 博客 | 论坛 | 百科 | GEC | 活动 | 主题月 | 电子展
返回列表 回复 发帖

一维小波变换,可多次分解

一维小波变换,可多次分解

1、题目:一维小波变换,可多次分解
2、原理:卷积核变为Daubechies正交小波基h[]和g[]的交替形式。增加了多次分解的功能。
3、代码:
  • #include <stdio.h>
  • #include <stdlib.h>
  • #include <math.h>
  • #define LENGTH 4096//信号长度
  • /*****************************************************************
  • * 一维卷积函数
  • *
  • * 说明: 循环卷积,卷积结果的长度与输入信号的长度相同
  • *
  • * 输入参数: data[],输入信号; h[],Daubechies小波基低通滤波器系数;
  • *            g[],Daubechies小波基高通滤波器系数; cov[],卷积结果;
  • *            n,输入信号长度; m,卷积核长度.
  • *
  • * 李承宇, lichengyu2345@126.com
  • *
  • *  2010-08-22
  • *****************************************************************/
  • void Covlution(double data[], double h[], double g[], double cov[]
  •                , int n, int m)
  • {
  •     int i = 0;
  •     int j = 0;
  •     int k = 0;
  •     //将cov[]清零
  •     for(i = 0; i < n; i++)
  •     {
  •         cov = 0;
  •     }
  •     //****************************************************
  •     //奇数行用h[]进行卷积
  •     //****************************************************
  •     //前m/2+1行
  •     i = 0;
  •     for(j = 0; j < m/2; j+=2, i+=2)
  •     {
  •         for(k = m/2-j; k < m; k++ )
  •         {
  •             cov += data[k-(m/2-j)] * h[k];//k针对core[k]
  •         }
  •         for(k = n-m/2+j; k < n; k++ )
  •         {
  •             cov += data[k] * h[k-(n-m/2+j)];//k针对data[k]
  •         }
  •     }
  •     //中间的n-m行
  •     for( ; i <= (n-m)+m/2; i+=2)
  •     {
  •         for( j = 0; j < m; j++)
  •         {
  •             cov += data[i-m/2+j] * h[j];
  •         }
  •     }
  •     //最后m/2-1行
  • //  i = ( (n - m) + m/2 + 1 )/2*2;//**********
  •     for(j = 1; j <= m/2; j+=2, i+=2)
  •     {
  •         for(k = 0; k < j; k++)
  •         {
  •             cov += data[k] * h[m-j-k];//k针对data[k]
  •         }
  •         for(k = 0; k < m-j; k++)
  •         {
  •             cov += h[k] * data[n-(m-j)+k];//k针对core[k]
  •         }
  •     }
  •     //****************************************************
  •     //偶数行用g[]进行卷积
  •     //****************************************************
  •     //前m/2+1行
  •     i = 1;
  •     for(j = 0; j < m/2; j+=2, i+=2)
  •     {
  •         for(k = m/2-j; k < m; k++ )
  •         {
  •             cov += data[k-(m/2-j)] * g[k];//k针对core[k]
  •         }
  •         for(k = n-m/2+j; k < n; k++ )
  •         {
  •             cov += data[k] * g[k-(n-m/2+j)];//k针对data[k]
  •         }
  •     }
  •     //中间的n-m行
  •     for( ; i <= (n-m)+m/2; i+=2)
  •     {
  •         for( j = 0; j < m; j++)
  •         {
  •             cov += data[i-m/2+j] * g[j];
  •         }
  •     }
  •     //最后m/2-1行
  • //  i = ( (n - m) + m/2 + 1 ) ;//*********
  •     for(j = 1; j <= m/2; j+=2, i+=2)
  •     {
  •         for(k = 0; k < j; k++)
  •         {
  •             cov += data[k] * g[m-j-k];//k针对data[k]
  •         }
  •         for(k = 0; k < m-j; k++)
  •         {
  •             cov += g[k] * data[n-(m-j)+k];//k针对core[k]
  •         }
  •     }
  • }
  • /*****************************************************************
  • *   排序函数
  • *
  • *   将卷积后的结果进行排序,使尺度系数和小波系数分开
  • *****************************************************************/
  • void Sort(double data[], double sort[], int n)
  • {
  •     for(int i = 0; i < n; i+=2)
  •     {
  •         sort[i/2] = data;
  •     }
  •     for(i = 1; i < n; i+=2)
  •     {
  •         sort[n/2+i/2] = data;
  •     }
  • }
  • /*****************************************************************
  • * 一维小波变换函数
  • *
  • * 说明: 一维小波变换,可进行多次分解
  • * 输入参数: input[],输入信号; output[],小波变换结果,包括尺度系数
  • * 和小波系数两部分; temp[],存放中间结果;h[],Daubechies小波基低通滤
  • * 波器系数;g[],Daubechies小波基高通滤波器系数;n,输入信号长度; m,
  • * Daubechies小波基紧支集长度; nStep,小波变换分解次数
  • *****************************************************************/
  • void DWT1D(double input[], double output[], double temp[], double h[],
  •            double g[], int n, int m, int nStep)
  • {
  •     int i = 0;
  •     for(i = 0; i < n; i++)
  •     {
  •         output = input;
  •     }
  •     for(i = 0; i < nStep; i++)
  •     {
  •         Covlution(output, h, g, temp, n, m);
  •         Sort(temp, output, n);
  •         n = n/2;
  •     fp=fopen("data.txt","r");
  •     if(fp==NULL) //如果读取失败
  •     {
  •         printf("错误!找不到要读取的文件/"data.txt/"/n");
  •         exit(1);//中止程序
  •     }
  •     while( fgets(s, 32, fp) != NULL )//读取长度n要设置得长一点,要保证读到回车符,这样指针才会定位到下一行?回车符返回的是零值?是,非数字字符经过atoi变换都应该返回零值
  •     {
  •     //  fscanf(fp,"%d", &data[count]);//一定要有"&"啊!!!最后读了个回车符!适应能力不如atoi啊
  •         data[n] = atof(s);
  •         n++;
  •     }
  •     //一维小波变换
  •     DWT1D(data, data_output, temp, h, g, n, m, nStep);
  •     //一维小波变换后的结果写入txt文件
  •     fp=fopen("test.txt","w");
  •     //打印一维小波变换后的结果
  •     for(i = 0; i < n/pow(2,nStep-1); i++)///pow(2,nStep-1)
  •     {
  •         printf("%f/n", data_output);
  •         fprintf(fp,"%f/n", data_output);
  •     }
  •     //关闭文件
  •     fclose(fp);
  • }

#include <stdio.h>#include <stdlib.h>#include <math.h>#define LENGTH 4096//信号长度/****************************************************************** 一维卷积函数** 说明: 循环卷积,卷积结果的长度与输入信号的长度相同** 输入参数: data[],输入信号; h[],Daubechies小波基低通滤波器系数;*  g[],Daubechies小波基高通滤波器系数; cov[],卷积结果;* n,输入信号长度; m,卷积核长度.** 李承宇, lichengyu2345@126.com** 2010-08-22 *****************************************************************/void Covlution(double data[], double h[], double g[], double cov[] , int n, int m){ int i = 0; int j = 0; int k = 0; //将cov[]清零 for(i = 0; i < n; i++) { cov = 0; } //**************************************************** //奇数行用h[]进行卷积 //**************************************************** //前m/2+1行 i = 0; for(j = 0; j < m/2; j+=2, i+=2) { for(k = m/2-j; k < m; k++ ) { cov += data[k-(m/2-j)] * h[k];//k针对core[k] } for(k = n-m/2+j; k < n; k++ ) { cov += data[k] * h[k-(n-m/2+j)];//k针对data[k] } } //中间的n-m行 for( ; i <= (n-m)+m/2; i+=2) { for( j = 0; j < m; j++) { cov += data[i-m/2+j] * h[j]; } } //最后m/2-1行// i = ( (n - m) + m/2 + 1 )/2*2;//********** for(j = 1; j <= m/2; j+=2, i+=2) { for(k = 0; k < j; k++) { cov += data[k] * h[m-j-k];//k针对data[k] } for(k = 0; k < m-j; k++) { cov += h[k] * data[n-(m-j)+k];//k针对core[k] } } //**************************************************** //偶数行用g[]进行卷积 //**************************************************** //前m/2+1行 i = 1; for(j = 0; j < m/2; j+=2, i+=2) { for(k = m/2-j; k < m; k++ ) { cov += data[k-(m/2-j)] * g[k];//k针对core[k] } for(k = n-m/2+j; k < n; k++ ) { cov += data[k] * g[k-(n-m/2+j)];//k针对data[k] } } //中间的n-m行 for( ; i <= (n-m)+m/2; i+=2) { for( j = 0; j < m; j++) { cov += data[i-m/2+j] * g[j]; } } //最后m/2-1行// i = ( (n - m) + m/2 + 1 ) ;//********* for(j = 1; j <= m/2; j+=2, i+=2) { for(k = 0; k < j; k++) { cov += data[k] * g[m-j-k];//k针对data[k] } for(k = 0; k < m-j; k++) { cov += g[k] * data[n-(m-j)+k];//k针对core[k] } }}/****************************************************************** 排序函数** 将卷积后的结果进行排序,使尺度系数和小波系数分开*****************************************************************/void Sort(double data[], double sort[], int n){ for(int i = 0; i < n; i+=2) { sort[i/2] = data; } for(i = 1; i < n; i+=2) { sort[n/2+i/2] = data; }}/****************************************************************** 一维小波变换函数** 说明: 一维小波变换,可进行多次分解 ** 输入参数: input[],输入信号; output[],小波变换结果,包括尺度系数* 和小波系数两部分; temp[],存放中间结果;h[],Daubechies小波基低通滤* 波器系数;g[],Daubechies小波基高通滤波器系数;n,输入信号长度; m,* Daubechies小波基紧支集长度; nStep,小波变换分解次数** 李承宇, lichengyu2345@126.com** 2010-08-22 *****************************************************************/void DWT1D(double input[], double output[], double temp[], double h[], double g[], int n, int m, int nStep){ int i = 0; for(i = 0; i < n; i++) { output = input; } for(i = 0; i < nStep; i++) { Covlution(output, h, g, temp, n, m); Sort(temp, output, n); n = n/2; }}void main(){ double data[LENGTH];//输入信号 double temp[LENGTH];//中间结果 double data_output[LENGTH];//一维小波变换后的结果 int n = 0;//输入信号长度 int m = 6;//Daubechies正交小波基长度 int nStep = 6;//分解级数 int i = 0; char s[32];//从txt文件中读取一行数据 static double h[] = {.332670552950, .806891509311, .459877502118, -.135011020010, -.085441273882, .035226291882}; static double g[] = {.035226291882, .085441273882, -.135011020010, -.459877502118, .806891509311, -.332670552950}; //读取输入信号 FILE *fp; fp=fopen("data.txt","r"); if(fp==NULL) //如果读取失败 { printf("错误!找不到要读取的文件/"data.txt/"/n"); exit(1);//中止程序 } while( fgets(s, 32, fp) != NULL )//读取长度n要设置得长一点,要保证读到回车符,这样指针才会定位到下一行?回车符返回的是零值?是,非数字字符经过atoi变换都应该返回零值 { // fscanf(fp,"%d", &data[count]);//一定要有"&"啊!!!最后读了个回车符!适应能力不如atoi啊 data[n] = atof(s); n++; } //一维小波变换 DWT1D(data, data_output, temp, h, g, n, m, nStep); //一维小波变换后的结果写入txt文件 fp=fopen("test.txt","w"); //打印一维小波变换后的结果 for(i = 0; i < n/pow(2,nStep-1); i++)///pow(2,nStep-1) { printf("%f/n", data_output); fprintf(fp,"%f/n", data_output); } //关闭文件 fclose(fp);}4、测试结果:
输入信号x(i)为:
取f1 = 5, f2 = 10, f0 = 320, n = 512。x(i)如图1所示:
图1 输入信号
各级分解的结果如图2~图7所示,左半部分为尺度系数,右半部分为小波系数:
图2 1级分解结果
图3 2级分解结果
图4 3级分解结果
图5 4级分解结果
图6 5级分解结果
图7 6级分解结果
图8是各级小波系数和第6级尺度系数的完整结果:
图8 第6级尺度系数和各级小波系数的完整结果
返回列表