本帖最后由 Nate_ 于 2016-4-17 15:57 編輯

points=1024 時，有波形輸出，但信號有5438個點。改為5438就不行。主程序：

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

close all;

points=5438;? ?? ???level=4;? ? sr=360;? ?num_inter=6;? ?wf='db4';

%所處理數據的長度 分解的級數? ?抽樣率? ???迭代次數? ?? ?小波名稱

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 [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

子程序2：

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)

interval(i+1)=intemin;

end

end

end

pc3inte=interval;

子程序3：

function??w2=Py_Pgama(w1,wpeak,wframe,level,sr)

% 該函數用于進行 Pgama 和 Py 投影

err=wpeak-w1.*(wpeak~=0);

w2=zeros(size(wpeak));

[r]=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

子程序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,-90000,90000]);

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

2016-4-17 15:52 上傳

2016-4-17 15:47 上傳

點擊文件名下載附件

70.54 KB, 下載次數: 232

信號

2016-4-17 15:48 上傳

點擊文件名下載附件

1.37 KB, 下載次數: 49

主程序

2016-4-17 15:48 上傳

點擊文件名下載附件

368 Bytes, 下載次數: 47

子程序

2016-4-17 15:48 上傳

點擊文件名下載附件

832 Bytes, 下載次數: 43

子程序

2016-4-17 15:49 上傳

點擊文件名下載附件

885 Bytes, 下載次數: 43

子程序

2016-4-17 15:49 上傳

點擊文件名下載附件

1.32 KB, 下載次數: 45

子程序