mjdzh.m

来自「这是一个用模极大值算法对一维信号去噪的程序,信号重构采用经典的交替投影法.」· M 代码 · 共 96 行

clear;close all;
load leleccum;
%[a1,a2]=textread('bc7.txt','%n%n','headerlines',2);
%x=a2(210:1233);
%x=a2(1:1024);
%x=x';

x=leleccum(1:1024);
sig=x;
x=awgn(x,10,30);
figure;
subplot(211);plot(sig);title('原始信号');grid on;
subplot(212);plot(x,'r');title('染噪信号');grid on;
points=length(x);level=4;sr=360;num_inter=6;wf='db4';
[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters('db4');
[swa,swd]=swt(x,4,Lo_D,Hi_D);
figure(2);
for i=1:level
    subplot(4,2,2*i-1);plot(swa(i,:));Ylabel(['近似系数A',num2str(i)]);
    subplot(4,2,2*i);plot(swd(i,:));Ylabel(['细节系数D',num2str(i)]);
end
%求小波变换的模极大值及其位置，并按层给出小波变换模极大值的波形
[ddw,wpeak]=hmjdz(swd,points,4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%进行模极大值的预处理
%设置阈值
c=1.8;
D4_wpeak=wpeak(level,:);
M=max(D4_wpeak);
thr=c*M/level;
%thr=10;%1/sqrt(2);
D4_wpeak=D4_wpeak.*(abs(D4_wpeak)>thr);
%在尺度上极大值点位置，构造一个搜索区域
%在尺度中，极大值点落在该区域的点保留，其他的置0
D3_wpeak=wpeak(level-1,:);
D4_p=(D4_wpeak~=0);
O_d4=2^(level-1);
for P_d4=O_d4:(length(D4_wpeak)-O_d4)
    if D4_p(P_d4)==1
        for i=1:O_d4-1
            D4_p(P_d4-i)=1;
        end
    end
end
D3_wpeak=D3_wpeak.*D4_p;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D2_wpeak=wpeak(level-2,:);
D3_p=(D3_wpeak~=0);
O_d3=2^(level-2);
for P_d3=O_d3:(length(D3_wpeak)-O_d3)
    if D3_p(P_d3)==1
        for i=1:O_d3-1
            D3_p(P_d3-i)=1;
        end
    end
end
D2_wpeak=D2_wpeak.*D3_p;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D1_wpeak=wpeak(1,:);
D2_p=(D2_wpeak~=0);
D1_wpeak=D1_wpeak.*D2_p;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wpeak=[D1_wpeak' D2_wpeak' D3_wpeak' D4_wpeak'];
wpeak=wpeak';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%重构信号
pswa=swa(level,:);
wframe=(wpeak~=0);
%迭代初始化
w0=zeros(1,points);
[a,d]=swt(w0,level,Lo_D,Hi_D);
w2=d;
for j=1:num_inter
    w2=Py_Pgama(d,wpeak,wframe,1,sr);
    w0=iswt(pswa,w2,Lo_R,Hi_R);
    [a,d]=swt(w0,level,Lo_D,Hi_D);
end
pswa=iswt(swa(level,:),w2,Lo_R,Hi_R);
figure;
subplot(311);plot(sig(1:points));title('原始信号');
subplot(312);plot(pswa(1:points));title('重构信号');
subplot(313);plot(x(1:points),'r');title('染噪信号');
werr=w2-swd;
figure;
for m=1:level
    wsnr(m)=20*log10(norm(swd(m,:))/norm(werr(m,:)));
    subplot(level+1,1,m);
    plot(w2(m,:),'r');grid on;Ylabel(strcat('j=',num2str(m)));
end
err=pswa(1:points)-x(1:points);
snr=20*log10(norm(x)/norm(err))