// 'haar' & 'morlet'
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "nrutil.h"
#define PI (3.141592653589793)
void wavefun(float *psi, float *xt, int wn, int len);
void morlet(float *psi, float *xt, float lb, float ub, float f0, int N);
void haar(float *psi, float *xt, int N);
void conv(float xi[], int nx, float yi[], int ny, float zo[]);
enum waveletname {
HAAR = 1, MORLET = 4 //保持其与matlab的相似性.
}wname;
float **cwt(float *sig, int m, float *scales, int n, char wname[])
// sig : 输入信号(长度为m)
// scales : 尺度向量(长度为n)
// wname : 小波母函数名字
{
const char sHAAR[] = "haar";
const char sMORLET[] = "morlet";
int precis = 10;
int i, j, xm, N, wn, len = pow(2, precis) + 2;
float xstep, xmax, a;
float *psi = vector(1, len);
float *xt = vector(1, len);
if (!strcmp(wname, sHAAR)) //strcmp函数返回字符差值,若两字符串相同则返回-0-.
wn = HAAR; //'haar'小波
else if (!strcmp(wname, sMORLET))
{
wn = MORLET; //'morlet'小波
len = len - 2;
}
else
{
wn = -1; //'-1'表示输入错误
return 0;
}
wavefun(psi, xt, wn, len);
xstep = xt[2] - xt[1];
xmax = xt[len]- xt[1];
float **coefs = matrix(1, n, 1, m);
for (i = 1; i <= n; i++)
{
a = scales[i];
N = floor(a * xmax);
/* #### 分配临时存储空间 #### */
float *tmp = vector(1, N+1);
float *wcoefs= vector(1, N+m);
/* ## 伸缩变换 ## */
for (j = 0; j <= N; j++)
{
xm = floor( j / (a * xstep) ) + 1;
tmp[ N+1 - j ] = psi[ xm ];
}
/* ## 卷积 ## */
conv(sig, m, tmp, N+1, wcoefs);
/* ## 微分 ## */
for (j = 1; j < N+m; j++)
wcoefs[j] = wcoefs[j+1] - wcoefs[j];
/* ## 数据截取 ## */
for (j = 1; j <= m; j++)
coefs[i][j] = wcoefs[ (N-1)/2 + j ] * -sqrt(a);
/* #### 释放临时存储空间 #### */
free_vector(tmp, 1, N+1);
free_vector(wcoefs, 1, N+m);
}
return coefs;
}
//#####################################################
void wavefun(float *psi, float *xt, int wn, int len)
// wn : wavelet name. len : length of psi.
{
int i;
float xstep;
switch (wn)
{
case HAAR:
// psi[1..len] = {0.0, 1.0, .. , 1.0, -1.0, .. , -1.0, 0};
haar(psi, xt, len);
break;
case MORLET:
// psi[1..len] = {..};
morlet(psi, xt, -8, 8, 5/(2*PI), len);
break;
default:
break;
}
// cumulative sum of psi.
xstep= xt[2] - xt[1];
for (i = 1; i <= len; i++)
{
psi[i] = psi[i] + psi[i-1];
psi[i-1] *= xstep;
}
psi[len] *= xstep;
}
// ############################################################################
// #### 小波母函数 ####
// ##1 'haar'小波
void haar(float *psi, float *xt, int N)
// psi : 小波母函数(时域)
// xt : 采样时刻点(t)
// N : 采样点数
{
int i;
xt[1] = psi[1] = 0.0;
for (i = 2; i <= N; i++)
{
psi[i] = (i > N/2) ? -1.0 : 1.0;
xt[i] = xt[i-1] + (double)1.0 / (N-2);
}
psi[N] = 0.0; //将psi[N]重新置零
}
// ##2 'morlet'小波 /* 实Morlet小波函数<时域>函数值 */
void morlet (float *psi, float *xt, float lb, float ub, float f0, int N)
// psi : 小波母函数(时域)
// xt : 采样时刻点(t)
// lb~ub : 取值范围(小波支撑集)
// f0 : 小波母函数的中心频率(2*PI*f0 >= 5)
// N : 采样点数
{
int i;
double siz = (ub - lb) / (N-1), t;
t = lb;
for(i = 1; i <= N; i++)
{
psi[i] = exp( -pow(t, 2) / 2 ) * cos(2 * PI * f0 * t); //实部
xt [i] = t;
t += siz;
}
}
// ############################################################################