[转帖]小波变换的模极大值重建信号源程序

AaronSpark

这是用小波变换模极大值重建信号的源程序，数据是一心电信号，在matlab6。5下实现，来源于胡广书的《现代信号处理教程》附属光盘，现提供给大家供大家学习参考，滤波部分可以根据个人情况进行修改。程序包含四个部份，其中wave_peak实现信号的分解并求出模极大值<BR><BR>%--------------------------------------------------------------------------<BR>% exa130301.m 例13.3.1: 利用小波变换模极大重建原信号<BR>%--------------------------------------------------------------------------<BR>close all;<BR>points=1024; level=6; sr=360; num_inter=6; wf='db3';<BR>%所处理数据的长度 分解的级数 抽样率 迭代次数 小波名称<BR><BR>offset=0;<BR>[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);<BR><BR>%计算小波分解系数和模极大序列<BR>[signal,swa,swd,ddw,wpeak]=wave_peak(points,level,Lo_D,Hi_D,Lo_R,Hi_R,offset);<BR>% signal: 原始信号; swa:小波概貌; swd:小波细节;<BR>% ddw: 局部极大位置; wpeak:小波变换的局部极大序列]<BR><BR>pswa=swa(level,: ); % pswa: 为待重建的信号<BR>wframe=(wpeak~=0);<BR>%迭代初始化<BR>w0=zeros(1,points);<BR>[a,d]=swt(w0,level,Lo_D,Hi_D);<BR>w2=d; % w2为待重建小波<BR>for j=1:num_inter<BR>w2=Py_Pgama(d,wpeak,wframe,1,sr); % 先进行Py投影和 Pgama投影<BR>w0=iswt(pswa,w2,Lo_R,Hi_R); % 再进行Pv投影<BR>[a,d]=swt(w0,level,Lo_D,Hi_D); % Pv<BR>end<BR>pswa=iswt(swa(level,: ),w2,Lo_R,Hi_R); % 计算重建信号<BR><BR>% 原信号和由模极大重建信号的比较<BR>figure,<BR>subplot(211)<BR>plot(pswa(1:points)); <BR>subplot(212)<BR>plot(signal(1:points),'r');<BR><BR>%分别计算重建小波以及原信号的信噪比<BR>werr=w2-swd;<BR>% 原信号的小波变换（swd)和重建后的小波变换（w2）的比较<BR>figure,<BR>for m=1:level<BR>wsnr(m)=20*log10(norm(swd(m,: ))/norm(werr(m,: )))<BR>subplot(level+1,1,m);<BR>plot(swd(m,: )),hold on,<BR>plot(w2(m,: ),'r');grid on;ylabel(strcat('j=',num2str(m))),axis tight;<BR>end<BR>err=pswa(1:points)-signal(1:points);<BR>snr=20*log10(norm(signal)/norm(err)) <BR><BR>

AaronSpark

function [signal,swa,swd,ddw,wpeak]=wave_peak(points,level,Lo_D,Hi_D,Lo_R,Hi_R,offset)<BR>% 该函数用于读取ecg信号，找到小波变换模极大序列<BR><BR>warning off;<BR>ecgdata=load('ecg.txt' );<BR>plot(ecgdata(1:points)),grid on,axis tight,axis([1,points,-50,300]);<BR>signal=ecgdata(1:points)'+offset;<BR><BR>% 信号的小波变换，按级给出概貌和细节的波形<BR>[swa,swd] = swt(signal,level,Lo_D,Hi_D);<BR>figure;<BR>subplot(level,1,1); plot(real(signal)); grid on;axis tight;<BR>for i=1:level<BR>subplot(level+1,2,2*(i)+1);<BR>plot(swa(i,: )); axis tight;grid on;xlabel('time');<BR>ylabel(strcat('a ',num2str(i)));<BR>subplot(level+1,2,2*(i)+2);<BR>plot(swd(i,: )); axis tight;grid on;<BR>ylabel(strcat('d ',num2str(i)));<BR>end<BR><BR>%求小波变换的模极大值及其位置<BR>ddw=zeros(size(swd));<BR>pddw=ddw;<BR>nddw=ddw;<BR>posw=swd.*(swd&gt;0);<BR>pdw=((posw(:,1:points-1)-posw(:,2:points))&lt;0);<BR>pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))&gt;0);<BR>negw=swd.*(swd&lt;0);<BR>ndw=((negw(:,1:points-1)-negw(:,2:points))&gt;0);<BR>nddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))&gt;0);<BR>ddw=pddw|nddw;<BR>ddw(:,1)=1;<BR>ddw(:,points)=1;<BR>wpeak=ddw.*swd;<BR>wpeak(:,1)=wpeak(:,1)+1e-10;<BR>wpeak(:,points)=wpeak(:,points)+1e-10;<BR><BR>%按级给出小波变换模极大的波形<BR>figure;<BR>for i=1:level<BR>subplot(level,1,i);<BR>plot(wpeak(i,: )); axis tight;grid on;<BR>ylabel(strcat('j= ',num2str(i)));<BR>end <BR>

sandy3lu

<P>thanks </P>

[em04]

cforlove

太好了,一直在找相关的东东

wangjiao

谢谢谢谢

q2384

bu不错

chinawang

想把这个程序做到LABVIEW里面去，但是不知道该怎样处理那两个子函数

谁能指点下    不胜感激

PSZ10

搂主你好，我现在用到二次样条小波变换及变换后极大，极小值的位置提取，可否介绍一下在那能找到源代码，请联系，谢谢。信箱：PSZ10@126.COM

pipi123

LZ没发原信号数据文件，我补上。

星光2006

谢谢提供,小波变换的模极大值重建信号源程序...

lsq001

这个例子很有代表性吧
支持一下

dreamseeker

呵呵  你发的这些代码还是很适用的，挺有代表性，要好好模仿下

johnbenny

但是这个例子有没有滤波的部分呢？

johnbenny

是不是通过模极大算法滤波的呢？

dewlily

好东西,支持LZ:@)