本代码是摘自胡广书编写的《现代信号处理教程》，对相关程序段进行分析说明。

主程序：

%小波模极大值重构是采用的交替投影法

close all;
points=1024;        level=6;    sr=360;   num_inter=6;   wf='db3';
%所处理数据的长度 分解的级数   抽样率     迭代次数      小波名称

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

%计算小波分解系数和模极大序列
[signal,swa,swd,ddw,wpeak]=wave_peak(points,level,Lo_D,Hi_D,Lo_R,Hi_R,offset);
% signal:  原始信号;       swa:小波概貌;  swd:小波细节;
% ddw:     局部极大位置; wpeak:小波变换的局部极大序列]

pswa=swa(level,:);  % pswa: 为待重建的信号
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);  % 先进行Py投影和 Pgama投影
       w0=iswt(pswa,w2,Lo_R,Hi_R);         % 再进行Pv投影
       [a,d]=swt(w0,level,Lo_D,Hi_D);      % Pv
    end
     pswa=iswt(swa(level,:),w2,Lo_R,Hi_R); % 计算重建信号
    
% 原信号和由模极大重建信号的比较
figure,
subplot(211)
plot(pswa(1:points));
subplot(212)
plot(signal(1:points),'r');

%分别计算重建小波以及原信号的信噪比
werr=w2-swd;
% 原信号的小波变换（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))

%用到的函数

（1）function pc3inte=P_y(interval,len);
% 该函数对区间进行裁减即Py投影，返回裁剪后的区间信号

if sign(interval(1))==sign(interval(len))
    interval=interval.*(sign(interval)==sign(interval(1)));
    inte=interp1([1,len],[interval(1),interval(len)],(1:len),'linear');
    interval=sign(interval(1))*(abs(inte)-(abs(inte)-abs(interval)).*((abs(inte)-abs(interval))>0));
else
    sgn=sign(interval(len)-interval(1));
    intemax=max([interval(1),interval(len)]);
    intemin=min([interval(1),interval(len)]);
    for i=1:len-2
        if sign(interval(i+1)-interval(i))~=sgn
            interval(i+1)=interval(i);
        end
        if interval(i+1)>intemax
            interval(i+1)=intemax;
        end
       if interval(i+1)<intemin
                interval(i+1)=intemin;
        end
    end
end
pc3inte=interval;

（2）function  w2=Py_Pgama(w1,wpeak,wframe,level,sr);
% 该函数用于进行 Pgama 和 Py 投影
err=wpeak-w1.*(wpeak~=0);
w2=zeros(size(wpeak));
[r,c]=size(wpeak);
% 对每一级小波分别进行处理
for m=1:r
    frame=find(wpeak(m,:));
    num_interval=length(frame)-1;
   
    % 先找到以模极大划分的区间, 然后对每一区间进行Py投影
      for j=1:num_interval
          interval=w1(m,frame(j):frame(j+1));
          len=length(interval);
             if len>2
                w1(m,frame(j):frame(j+1))=P_y(interval,len);
             end
      end
  
      % 再逐一区间进行Pgama投影
      for j=1:num_interval
          interval=err(m,frame(j):frame(j+1));
          if r==1
              err(m,frame(j):frame(j+1))=P_gama(interval,level,sr);
          else
              err(m,frame(j):frame(j+1))=P_gama(interval,m,sr);
          end
      end
      w2(m,:)=w1(m,:)+err(m,:);
end

（3）function [inter]=P_gama(interval,lev,sr);
T=length(interval);
%该函数对一个区间进行Pgama投影，返回修正的区间
if T==2
    inter=interval;
else
    t=linspace(0,(T-1)/sr,T);
    para=(([1,1;exp(2^(-lev)*t(T)),exp(-2^(-lev)*t(T))])\[interval(1),interval(T)]')';
    alpha=para(1);
    beta=para(2);
    inter=alpha.*exp(2^(-lev).*t)+beta.*exp(-2^(-lev).*t);
end

（4）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 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;

%按级给出小波变换模极大的波形
figure;
for i=1:level
    subplot(level,1,i);
    plot(wpeak(i,:)); axis tight;grid on;
ylabel(strcat('j=   ',num2str(i)));
end

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