图像处理之基础---一维小波变换，可多次分解

时间：2014-09-01 00:14:12 阅读：361 评论：0 收藏：0 [点我收藏+]

1、题目：一维小波变换，可多次分解

2、原理：卷积核变为Daubechies正交小波基h[]和g[]的交替形式。增加了多次分解的功能。

3、代码：

[cpp] view plain copy

#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[i] = 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[i] += data[k-(m/2-j)] * h[k];//k针对core[k]
}
for(k = n-m/2+j; k < n; k++ )
{
cov[i] += 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[i] += 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[i] += data[k] * h[m-j-k];//k针对data[k]
}
for(k = 0; k < m-j; k++)
{
cov[i] += 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[i] += data[k-(m/2-j)] * g[k];//k针对core[k]
}
for(k = n-m/2+j; k < n; k++ )
{
cov[i] += 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[i] += 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[i] += data[k] * g[m-j-k];//k针对data[k]
}
for(k = 0; k < m-j; k++)
{
cov[i] += 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[i];
}
for(i = 1; i < n; i+=2)
{
sort[n/2+i/2] = data[i];
}
}
/*****************************************************************
*　一维小波变换函数
*
*　说明: 一维小波变换,可进行多次分解　
*
*　输入参数: 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[i] = input[i];
}
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[i]);
fprintf(fp,"%f/n", data_output[i]);
}
//关闭文件
fclose(fp);
}