%--------------------------------------------------------------------------
% exa130301.m  例13.3.1: 利用小波變換模極大重建原信號
%--------------------------------------------------------------------------
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))