這是用小波變換模極大值重建信號(hào)的源程序，數(shù)據(jù)是一心電信號(hào)，在matlab6。5下實(shí)現(xiàn)，來(lái)源于胡廣書(shū)的《現(xiàn)代信號(hào)處理教程》附屬光盤(pán)，現(xiàn)提供給大家供大家學(xué)習(xí)參考，濾波部分可以根據(jù)個(gè)人情況進(jìn)行修改。程序包含四個(gè)部份，其中wave_peak實(shí)現(xiàn)信號(hào)的分解并求出模極大值

%--------------------------------------------------------------------------
% exa130301.m 例13.3.1: 利用小波變換模極大重建原信號(hào)
%--------------------------------------------------------------------------
close all;
points=1024; level=6; sr=360; num_inter=6; wf='db3';
%所處理數(shù)據(jù)的長(zhǎng)度 分解的級(jí)數(shù) 抽樣率 迭代次數(shù) 小波名稱(chēng)

offset=0;
[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);

%計(jì)算小波分解系數(shù)和模極大序列
[signal,swa,swd,ddw,wpeak]=wave_peak(points,level,Lo_D,Hi_D,Lo_R,Hi_R,offset);
% signal: 原始信號(hào); swa:小波概貌; swd:小波細(xì)節(jié);
% ddw: 局部極大位置; wpeak:小波變換的局部極大序列]

pswa=swa(level,: ); % pswa: 為待重建的信號(hào)
wframe=(wpeak~=0);
%迭代初始化
w0=zeros(1,points);
[a,d]=swt(w0,level,Lo_D,Hi_D);
w2=d; % w2為待重建小波
for j=1:num_inter
w2=Py_Pgama(d,wpeak,wframe,1,sr); % 先進(jìn)行Py投影和 Pgama投影
w0=iswt(pswa,w2,Lo_R,Hi_R); % 再進(jìn)行Pv投影
[a,d]=swt(w0,level,Lo_D,Hi_D); % Pv
end
pswa=iswt(swa(level,: ),w2,Lo_R,Hi_R); % 計(jì)算重建信號(hào)

% 原信號(hào)和由模極大重建信號(hào)的比較
figure,
subplot(211)
plot(pswa(1:points)); 
subplot(212)
plot(signal(1:points),'r');

%分別計(jì)算重建小波以及原信號(hào)的信噪比
werr=w2-swd;
% 原信號(hào)的小波變換（swd)和重建后的小波變換（w2）的比較
figure,
for m=1:level
wsnr(m)=20*log10(norm(swd(m,: ))/norm(werr(m,: )))
subplot(level+1,1,m);
plot(swd(m,: )),hold on,
plot(w2(m,: ),'r');grid on;ylabel(strcat('j=',num2str(m))),axis tight;
end
err=pswa(1:points)-signal(1:points);
snr=20*log10(norm(signal)/norm(err))
function [signal,swa,swd,ddw,wpeak]=wave_peak(points,level,Lo_D,Hi_D,Lo_R,Hi_R,offset)
% 該函數(shù)用于讀取ecg信號(hào)，找到小波變換模極大序列

warning off;
ecgdata=load('ecg.txt' );
plot(ecgdata(1:points)),grid on,axis tight,axis([1,points,-50,300]);
signal=ecgdata(1:points)'+offset;

% 信號(hào)的小波變換，按級(jí)給出概貌和細(xì)節(jié)的波形
[swa,swd] = swt(signal,level,Lo_D,Hi_D);
figure;
subplot(level,1,1); plot(real(signal)); grid on;axis tight;
for i=1:level
subplot(level+1,2,2*(i)+1);
plot(swa(i,: )); axis tight;grid on;xlabel('time');
ylabel(strcat('a ',num2str(i)));
subplot(level+1,2,2*(i)+2);
plot(swd(i,: )); axis tight;grid 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;

%按級(jí)給出小波變換模極大的波形
figure;
for i=1:level
subplot(level,1,i);
plot(wpeak(i,: )); axis tight;grid on;
ylabel(strcat('j= ',num2str(i)));
end