C语言怎么实现小波变换

#include <stdio.h>

#include <stdlib.h>

#define LENGTH 512//信号长度

/******************************************************************

*　一维卷积函数

*

*　说明: 循环卷积,卷积结果的长度与输入信号的长度相同

*

*　输入参数: data[],输入信号 core[],卷积核 cov[],卷积结果

*  n,输入信号长度 m,卷积核长度.

*

*　李承宇, lichengyu2345@126.com

*

*  2010-08-18  

******************************************************************/

void Covlution(double data[], double core[], 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

}

//前m/2+1行

i = 0

for(j = 0 j < m/2 j++, i++)

{

for(k = m/2-j k < m k++ )

{

cov[i] += data[k-(m/2-j)] * core[k]//k针对core[k]

}

for(k = n-m/2+j k < n k++ )

{

cov[i] += data[k] * core[k-(n-m/2+j)]//k针对data[k]

}

}

//中间的n-m行

for( i = m/2 i <= (n-m)+m/2 i++)

{

for( j = 0 j < m j++)

{

cov[i] += data[i-m/2+j] * core[j]

}

}

//最后m/2-1行

i = (n - m) + m/2 + 1

for(j = 1 j < m/2 j++, i++)

{

for(k = 0 k 凯庆碧< j k++)

{

cov[i] += data[k] * core[m-j-k]//k针对data[k]

}

for(k = 0 k < m-j k++)

{

cov[i] += core[k] * data[n-(m-j)+k]//k针对core[k]

}

}

}

/******************************************************************

*　一维小波变换函数

*

*　说明: 一维小波变换,只变换一次

*

*　输入参数: input[],输入信号 output[],小波变换结果，包括尺度系数和

*　小波系数两部分 temp[],存放中间结果h[],Daubechies小波基低通滤波器系数

*　g[],Daubechies小波基高通滤波器系数n,输入信号长度 m,Daubechies小波基紧支集长度.

*

*　李承宇, lichengyu2345@126.com

*

*  2010-08-19  

******************************************************************/

void DWT1D(double input[], double output[], double temp[], double h[], 

   double g[], int n, int m)

{

// double temp[LENGTH] = {0}//?????????????

int i = 0

/*

//尺度系数和小波系数放在一起

Covlution(input, h, temp, n, m)

for(i = 盯举0 i < n i += 2)

{

output[i] = temp[i]

}

Covlution(input, g, temp, n, m)

for(i = 1 i < n i += 差烂2)

{

output[i] = temp[i]

}

*/

//尺度系数和小波系数分开

Covlution(input, h, temp, n, m)

for(i = 0 i < n i += 2)

{

output[i/2] = temp[i]//尺度系数

}

Covlution(input, g, temp, n, m)

for(i = 1 i < n i += 2)

{

output[n/2+i/2] = temp[i]//小波系数

}

}

void main()

{

double data[LENGTH]//输入信号

double temp[LENGTH]//中间结果

double data_output[LENGTH]//一维小波变换后的结果

int n = 0//输入信号长度

int m = 6//Daubechies正交小波基长度

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)

//一维小波变换后的结果写入txt文件

fp=fopen("data_output.txt","w")

//打印一维小波变换后的结果

for(i = 0 i < n i++)

{

printf("%f/n", data_output[i])

fprintf(fp,"%f/n", data_output[i])

}

//关闭文件

fclose(fp)

}

(1)小波模极大值重构 MATLAB代码_天天向上_新浪博客 http://blog.sina.com.cn/s/blog_6c00b0e30100u1s5.html

function

[signal,swa,swd,ddw,wpeak]=wave_peak(points,level,Lo_D,Hi_D,Lo_R,Hi_R,offset)

%

该函数用于读取ecg信号，派升找到小波变换模极大序列

warning off

ecgdata=load('ecg.txt')

%需要分析的信号，自己加

plot(ecgdata(1:points)),grid on,axis

tight,axis([1,points,-50,300])

signal=ecgdata(1:points)'+offset

% 信号的小波变换，按级给出概貌和细节的波形

[swa,swd] =

swt(signal,level,Lo_D,Hi_D)

figure

subplot(level,1,1)

plot(real(signal))grid onaxis tight

for i=1:level

subplot(level+1,2,2*(i)+1)

plot(swa(i,:))axis

tightgrid onxlabel('time')

ylabel(strcat('a

',num2str(i)))

subplot(level+1,2,2*(i)+2)

plot(swd(i,:))axis

tightgrid on

ylabel(strcat('d ',num2str(i)))

end

%求小波变换的模极大值及其位置

ddw=zeros(size(swd))

pddw=ddw

nddw=ddw

posw=swd.*(swd>0)

pdw=((posw(:,1:points-1)-posw(:,2:points))<0)

pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0)

negw=swd.*(swd<0)

ndw=((negw(:,1:points-1)-negw(:,2:points))>铅羡橡0)

nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0)

ddw=pddw|nddw

ddw(:,1)=1

ddw(:,points)=1

wpeak=ddw.*swd

wpeak(:,1)=wpeak(:,1)+1e-10

wpeak(:,points)=wpeak(:,points)+1e-10

%按级给出小波变换模极槐旁大的波形

figure

for i=1:level

subplot(level,1,i)

plot(wpeak(i,:))axis tightgrid

on

ylabel(strcat('j= ',num2str(i)))

end

注：运行此程序时一定要将待处理信号添加进去，程序中的红色部分。

追问：

ecgdata=load('ecg.txt')

---