本帖最后由 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
子程序
