#include <stdlib.h>
#define LENGTH 512//信号长度
/******************************************************************
* 一维卷积函数
*
* 说明: 循环卷积,卷积结果的长度与输入信号的长度相同
*
* 输入参数: data[],输入信号 core[],卷积核 cov[],卷积结果
* n,输入信号长度 m,卷积核长度.
*
* 李承宇, [email protected]
*
* 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小波基紧支集长度.
*
* 李承宇, [email protected]
*
* 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)
}
我只听闻小波技术，从未试过哈，只做过傅里叶变换而已。。。但是，你可以分析下你输出的数据，正确与否，就用matlab的输出和C执行器的输出对比下，就知道靠不靠谱了啊。如果嫌数据太多，可以找一些对比分析软件，如果实在无法全部分析，就弄几个特殊点，让程序就在那一些点输出对应值，对比matlab和C执行器的输出，多次换点，再对比，就知道靠不靠谱了。。。。你只是通过运算时间，怎么知道靠不靠谱呢，呵呵，要是你用多核DSP来算，只是一瞬间的事情，岂不是更不靠谱了。
(1)小波模极大值重构 MATLAB代码_天天向上_新浪博客 http://blog.sina.com.cn/s/blog_6c00b0e30100u1s5.htmlfunction
[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')
