小波变换模极大值去噪

小波变换模极大值去噪Matlab程序源码 包含模极大值提取函数、Pgama、Pv投影函数。

 N=1024;n1=7;
   n=0:N-1;
   Fs=500;
   t=n/Fs;
   x=1*sin（2*pi*10*t）;
   y=1.4*rand（11024）;
   for  i=1:length（y）
   if y（1i）>0.96
       y（i）=4;
       else if y（1i）>0.92
               y（1i）=2;
        else
               y（1i）=0;
   end
   end
   end
   a=0;
   for  i=1:length（y）
   if y（1i）>0
       a=a+1;
   else a=a;
   end
   end
   a
  signal=x+y;
  xref=x;
  points=1024;        level=4;    sr=360;   num_inter=6;   wf=‘db3‘; 
%所处理数据的长度    分解的级数   抽样率    迭代次数        小波名称 
offset=0; 
%____进行二进制小波变换（离散平稳小波变换），并给出各级波形： 
[Lo_DHi_DLo_RHi_R]=wfilters（wf）; 
[swaswd] = swt（signallevelLo_DHi_D）; 
figure; 
subplot（level11）; plot（real（signal））;
title（‘含噪声信号‘）grid on;axis tight; 
for i=1:level 
    subplot（level+122*（i）+1）; 
    plot（swa（i:））; axis tight;grid on;xlabel（‘time‘）; 
    ylabel（strcat（‘a   ‘num2str（i）））; 
    subplot（level+122*（i）+2）; 
    plot（swd（i:））; axis tight;grid on; 
ylabel（strcat（‘d   ‘num2str（i）））; 
end 
 
 
%____求小波变换的模极大值及其位置并按级给出小波变换模极大的波形: 
% swa:小波概貌;  swd:小波细节; 
% ddw:局部极大位置; wpeak:小波变换的局部极大序列。 
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; 
subplot（level+111）; plot（real（signal））; grid on;axis tight; 
for i=1:level 
    subplot（level+11i+1）; 
    plot（wpeak（i:））; axis tight;grid on; 
ylabel（strcat（‘j=   ‘num2str（i）））; 
end 
 
 
%____进行模极大值的处理： 
C=0.8;  
%此参数需要调节，为了在最大尺度上设定合适阈值，以确定最大尺度上该保留的模极大值点。 
D4_wpeak=wpeak（level:）; 
M=max（D4_wpeak）; 
Thr=C*M/level; %阈值计算，可参考论文：“3mm波段脉冲雷达系统研究和小波去噪分析“。 
D4_wpeak=D4_wpeak.*（abs（D4_wpeak）>Thr）; 
 
%模极大值的处理方式： 
%在尺度j上极大值点位置，构造一个搜索区域， 
%在尺度j-1中，极大值点落在该区域的点保留，其他的置0； 
D3_wpeak=wpeak（level-1:）; 
D4_p=（D4_wpeak~=0）; 
O_d4=3;%该参数确定在上一级搜索极大值的范围，可以调整。 
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=3;%该参数确定在上一级搜索极大值的范围，可以调整。 
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:）; % pswa: 为待重建的信号 
wframe=（wpeak~=0）; %迭代初始化 
w0=zeros（1points）; 
[ad]=swt（w0levelLo_DHi_D）; 
w2=d;  % w2为待重建小波 
    for j=1:num_inter 
       w2=Py_Pgama（dwpeakw

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       5025  2018-03-05 10:01  Module_maximum\maichongsin2.m

     文件        932  2018-03-05 10:01  Module_maximum\Py_Pgama.m

     文件        369  2003-12-03 17:09  Module_maximum\P_gama.m

     文件        847  2018-03-05 10:01  Module_maximum\P_y.m

     目录          0  2018-03-05 10:05  Module_maximum

----------- ---------  ---------- -----  ----

                 7173                    5