设为首页收藏本站|繁體中文 快速切换版块

 找回密码
 立即加入
搜索
查看: 3691|回复: 5

为什么模极大值法无法检测奇异点

[复制链接]

该用户从未签到

尚未签到

发表于 2009-4-23 08:13:56 | 显示全部楼层 |阅读模式

马上加入,结交更多好友,共享更多资料,让你轻松玩转电力研学社区!

您需要 登录 才可以下载或查看,没有账号?立即加入

×
原始图形如下,添加了噪声,使用的是db3小波,做了5层分解。
+ d/ w0 Z9 P+ r, r- g可是使用模极大值法检测不到第二个波的波头,不知道为什么?4 k8 n. y0 H* T" c; O1 D7 o+ ^
请达人指教。
: h: n) _: a9 S& B% x2 S' j源程序在如下% s8 H: t$ Y0 Y9 W: p
plot(cable(:,1),cable(:,2))0 {& ^- k0 }1 |$ ?6 i; y- t
tcable=cable(:,1);
3 w' n1 _. {( u# Dycable=cable(:,2);
# ?5 Z- |( B* w* k* [length(ycable)$ I0 n. y# R4 N% r/ k( a

+ f/ X; Z8 Q  f" l& p' |/ rsignal_length=33856;7 z! A5 m. _7 k1 I* |' R
t(1:signal_length,1)=0;
' Y/ x" u3 A6 F: a Ncable=length(tcable);" D4 {% p" y- l' P  i
for i=1:Ncable
' C5 a$ e- U6 k+ x( @4 L: L. G t(i,1)=tcable(i,1);
5 t3 o3 x; a' @  e- s6 S/ U/ I. p end# j" N. |) f2 u& T4 a, E2 X

# x- J; a* m$ K3 W1 b4 G' e- B& y& Vy(1:signal_length,1)=0;
$ v' p2 b+ l# @6 k6 Vfor i=1:Ncable' y/ o) k$ W) [7 q/ f0 S0 `
y(i,1)=ycable(i,1);
/ H7 J! E) }: \/ uend
8 B- m; s9 S% h noise=0.01*rand(1,length(y));) M4 @% }( j( C' c! `5 u
y_noise=y+noise';3 J! \! I- N! M" {: V# r, C0 ~
) \5 @! x- a. E6 P. q6 K
signal=y_noise;
; _$ s+ U  _6 b%plot(t,signal)% r! [& P6 _0 K* G& l; a
points=length(signal);        level=5;    sr=360;   num_inter=6;   wf='db3';! g( y7 r* P) L* `' W7 }! L
%所处理数据的长度    分解的级数   抽样率    迭代次数        小波名称( K% H- w, {$ v; ^
offset=0;8 p' K# y9 U) F0 V7 y: b& B- J
! d1 K6 ^( _! y! E3 \  W
! N( L. [7 D" c7 ~
%____进行二进制小波变换(离散平稳小波变换),并给出各级波形:
$ l' q5 K& I/ o. J) E& Z[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters(wf);%Lo_D:分解低通滤波器$ l9 V& `# M# _
                                   %Hi_D:分解高通滤波器
& e* k7 r' L; J6 i0 L+ S% C                                   %Lo_R:重建低通滤波器! n3 D+ c! {  e* w6 D/ {, G: o) t
                                   %Hi_R:重建高通滤波器( ]7 ]# n/ @# }& p0 G, L  E( W
                                   " O0 w, D6 a* C! e  O% b
[swa,swd] = swt(signal,level,Lo_D,Hi_D);%swa小波概貌; Z& }% P4 i* d1 H4 C
                                        %swd小波系数
) g: q% |6 C$ b) Z% ]* \5 S- @ " M# d) G9 I2 x) @
figure;%figure1
( S8 i: K' h& Q# c  }2 h0 ssubplot(level,1,1); plot(real(signal)); grid on;axis tight;
3 ~$ m! k: Y5 G% ~% V: a" O& E( E3 k: }for i=1:level  F. }9 c' f7 C0 {5 w
    subplot(level+1,2,2*(i)+1);# @5 W7 Z# U  p/ I$ I2 z0 M
    plot(swa(i,:)); axis tight;grid on;xlabel('time');  z. q2 Q9 B' F# y) s
    ylabel(strcat('a   ',num2str(i)));
: v4 d0 S' Z4 v1 ]- s' f    subplot(level+1,2,2*(i)+2);7 P: v9 S% T8 ]* a3 ^0 j4 V( P% Q: [
    plot(swd(i,:)); axis tight;grid on;
' A6 D( j) S* b; P' O' D" x1 Aylabel(strcat('d   ',num2str(i)));
" a! R6 h! N. u( B& O2 D8 Cend! `) I  n" R8 K! O' }
%以上内容是对信号进行离散小波变换。! |: T0 _# N& q6 h
%尺度为5
5 m4 C/ W2 z# B& k' D$ C& l! y, E& l( q4 I4 p4 j* w. [
- I* f5 O) k1 x/ u. N
8 M9 W, w( I+ ~
%____求小波变换的模极大值及其位置,并按级给出小波变换模极大的波形:9 A! G+ M, b% R/ t0 A
% swa:小波概貌;  swd:小波细节;; y3 `8 {% r$ ?' w" b0 S; @# P
% ddw:局部极大位置; wpeak:小波变换的局部极大序列。  G2 m0 }/ T- P. g
ddw=zeros(size(swd));%建立与swd维数相同的矩阵,矩阵值都为零
$ l& A. R* W1 k# i/ l8 V5 npddw=ddw;
9 ^$ u$ |8 A! n( d# Inddw=ddw;. p) D& p8 b" C6 h1 B
posw=swd.*(swd>0);%.*数组乘法,将swd中小于零的数置零
, @. R# y' A2 f  }0 u0 G6 W- B8 Q& lpdw=((posw(:,1:points-1)-posw(:,2:points))<0);" F/ _4 H$ K# k* L" ~9 z
pddw(:,2:points-1)=((pdw(:,1:points-2)-pdw(:,2:points-1))>0);) H2 ]  }* e3 F7 X1 P
negw=swd.*(swd<0);%将swd中大于零的数置零
8 e/ K* ~$ }  y0 ^- X0 Endw=((negw(:,1:points-1)-negw(:,2:points))>0);
6 r! V" k3 {+ bnddw(:,2:points-1)=((ndw(:,1:points-2)-ndw(:,2:points-1))>0);
6 ^8 t  G, w+ T1 |) Hddw=pddw|nddw;%|:或
$ v' c! x  U$ C  P) I- Jddw(:,1)=1;
4 ], Y) F2 |/ Y+ z6 lddw(:,points)=1;3 u- j" m; E# V# Z$ i' E5 U2 f1 h/ B$ Y
wpeak=ddw.*swd;0 {2 M4 ~! O7 w5 U& Q5 x0 k
wpeak(:,1)=wpeak(:,1)+1e-10;%第一列的值加上1e-10,不知道出于什么目的,目前看来不会对结果产生影响
' F: ], r3 {6 |. [wpeak(:,points)=wpeak(:,points)+1e-10;%同理,将最后一列的值加上1e-10,不知道出于什么目的,目前看来同样不会对结果产生影响* T7 S, T8 p2 P6 ~2 y7 `/ E* r3 A
%以上内容为:寻找每级尺度上小波变换系数对应的模极大值点6 c( e. G9 S( H! m3 ]$ c. ^' [  ^

, R  X& r2 e' N) K. _wpeak_threshold=wpeak(level,:)/max(abs(wpeak(level,:)));* {( C; s$ G! p$ ~2 J) Q' G
threshold=0.1;  %  阈值 & k. M& V  n/ V/ G3 Y* Q
, R0 |. ^/ e- L$ j& q, _( m
for i=1:points;0 W( M# m' @' E  K6 \; m  K$ p  }2 M
   if (abs(wpeak_threshold(i))<threshold)! B& Y7 j. H  n4 a  ~/ i& ]
        wpeak(level,i)=0;
9 e9 X9 d* P3 b' k, H2 o. F' A    end
7 e0 z& U7 y4 |1 T4 s& e: Dend;
  J- E; E8 X: L9 f9 a# c* P$ U+ K4 @! xD5_wpeak=wpeak(level,:);! x0 k) `2 }/ |5 K% s
3 S) a; e) O: d4 n1 N0 F
%对最大尺度2.^J上的极大值点设定阈值T0,将低于阈值T0的模极大值点去掉,3 s% ~3 s. v3 ^9 H2 t1 `

+ O  R, x8 |0 f
0 M4 x9 f* ]4 F  T5 {  {  I%____进行模极大值的处理:
& b( ~' e9 b9 Q8 p+ e3 U7 A%%C=0.8; # ?. j8 [& o+ k$ Q! P
%此参数需要调节,为了在最大尺度上设定合适阈值,以确定最大尺度上该保留的模极大值点。
# n: t6 ^0 z! ~' R, @' [0 A% q%%D5_wpeak=wpeak(level,:);4 j7 L# |- a0 F% [& c* S
%%M=max(abs(D5_wpeak));2 V! g( B% w( s3 `7 A
%%Thr=C*M/level; %阈值计算,可参考论文:"3mm波段脉冲雷达系统研究和小波去噪分析"。
  C$ e/ `  j1 j& X! d+ ]1 U6 m%%for i=1:points
) f4 W7 B; R) s" }& W8 m4 [%%    if(abs(D5_wpeak(i))<Thr)
( I1 @( W  I+ z; L. S9 S6 R8 Q' |%%        wpeak(level,i)=0;, o4 [  u: x9 [7 L2 l/ Z  g$ D( M
%%    end
8 e5 X" H& D+ C& _" d%%end
% ?# |3 J) r5 H. i4 P; [7 H, Z%%D5_wpeak=wpeak(level,:);
  x2 u5 P) O0 p2 `4 z4 W5 A: `& a! z- [: _9 R8 L
figure;%figure2  g* v8 ~2 V: m' F  _, J( l
subplot(level+1,1,1); plot(real(signal)); grid on;axis tight;
, Z0 g7 o( Z) _+ Ffor i=1:level
* X1 m+ M" R$ I( E2 H# A+ k3 T( ^    subplot(level+1,1,i+1);- V2 _7 u! m3 w0 I+ @
    plot(wpeak(i,:)); axis tight;grid on;
4 z3 A4 I4 }8 v5 t* ]9 Wylabel(strcat('j=   ',num2str(i)));
9 G2 {) [2 p% g8 Uend
* h/ i+ h$ X& P% M$ R, I/ I! Y$ Z! b& V, v( S5 Z
6 e. Z. \8 e4 e. [: P
%步骤4
7 }4 v5 j" u4 m- ]7 w0 D) t( M; g%模极大值的处理方式:
& }/ u- t& v. U, r* O3 S) K%在尺度j上极大值点位置,构造一个搜索区域,7 x! T5 {+ h' `4 Y7 E" F
%在尺度j-1中,极大值点落在该区域的点保留,其他的置0;
( k: b2 Z  w+ F: ?1 ^3 _D4_wpeak=wpeak(level-1,:);
( l2 z: Z+ r& \' f+ v% Qsousu=6;
+ Y  b' [4 M; V4 J6 ^0 ]; wD5_p=(D5_wpeak~=0);
7 S( F- A+ D8 U$ J/ |O_d5=sousu;%该参数确定在上一级搜索极大值的范围,可以调整。
8 K- ]: M! F& b! @3 ]' K1 U. P3 Pfor P_d5=O_d5:(length(D5_wpeak)-O_d5);8 g( `( N/ k; f; d  h9 I0 w
    if D5_p(P_d5)==1;
/ @2 ], V! E. S; n        for i=1:O_d5-1;
! \% j2 E+ S8 z! m        D5_p(P_d5-i)=1;  8 r+ k. S3 V% p/ O
        end ;6 M$ P. X, g2 b8 {$ `  U8 V
    end;     
9 f' g" z! r- b' r5 g9 f  q  F3 Q( Jend;
! N8 z/ `& C# f' v% K/ kD4_wpeak=D4_wpeak.*D5_p;
2 ]! O( D, G9 J) L6 F4 P, E, c. ]1 d5 W# N1 X: e
D3_wpeak=wpeak(level-2,:);* u( O* P! m" L: O2 a
D4_p=(D4_wpeak~=0);
, l+ G- \3 |* y9 MO_d4=sousu;%该参数确定在上一级搜索极大值的范围,可以调整。
3 s  p; E7 X+ M6 t/ D0 Q: Efor P_d4=O_d4:(length(D4_wpeak)-O_d4);( t/ x) ?+ b6 C/ n, r6 w, w. i8 g; Y5 ^
    if D4_p(P_d4)==1; / y2 Q9 X, T: q' [) ]1 m
        for i=1:O_d4-1;
8 _5 d: A# Z& o0 G& y        D4_p(P_d4-i)=1;  
* i5 O9 P/ E1 w  y3 q8 H        end ;
* P) p; ]7 f2 ]* j    end;     
! w3 h" F% v! A# S6 K! Hend;
6 W1 F$ G: E* z" nD3_wpeak=D3_wpeak.*D4_p;& ?5 e" N" Q9 I' t

, n+ U! [5 Z9 [& _D2_wpeak=wpeak(level-3,:);
( m/ a% a. o! q! wD3_p=(D3_wpeak~=0);0 N/ M3 I' B- K6 g
O_d3=sousu;%该参数确定在上一级搜索极大值的范围,可以调整。
/ ?% z5 I' N: b. Q' v* sfor P_d3=O_d3:(length(D3_wpeak)-O_d3);
: e  B! _4 w2 ~3 w$ ~8 N+ L    if D3_p(P_d3)==1;
' f; Z. A+ W* }1 t* {3 O0 W        for i=1:O_d3-1;* @" N$ r5 P$ Z, J
        D3_p(P_d3-i)=1;' n6 Y1 r$ `4 Z0 l8 Z8 t
        end ;. G# J0 K# b, m3 w" k7 b
    end;     + s$ _) X- T7 q5 N: J/ |
end;
: E3 H7 k$ M9 AD2_wpeak=D2_wpeak.*D3_p;
( X) W1 n7 X; ~: o1 \5 ~6 T, R% H' D
* A2 A" @; E8 s$ l, D' R/ r, K%第一层单独处理,在第二层极大值点位置上,保留第一层相应极大值点:+ M2 v. t# Q/ d% L
D1_wpeak=wpeak(level-4,:);
7 B+ \, {; L6 P. K. I8 RD2_p=(D2_wpeak~=0);9 Z9 B4 V+ `4 g3 \# n
D1_wpeak=D1_wpeak.*D2_p;, c3 K, X$ n! `& K

! h. h$ Z3 [" {1 q# Iwpeak=[D1_wpeak' D2_wpeak' D3_wpeak' D4_wpeak' D5_wpeak'];
" V" H0 e% W: P4 O; [6 Cwpeak=wpeak';
/ b  W# A* o: P9 t/ H5 ?
& b4 N  O. R9 T2 c7 u, t%____重构信号:
% L9 N0 h8 Z) G& Epswa=swa(level,:); % pswa: 为待重建的信号
3 m, S+ O  V8 I- hwframe=(wpeak~=0); %迭代初始化
* w% |& g, P6 m, I# ww0=zeros(1,points);+ s1 i5 w* i: T1 f# j2 D
[a,d]=swt(w0,level,Lo_D,Hi_D);& M( s$ K! E$ o" m$ i$ l& H/ G
w2=d;  % w2为待重建小波9 W+ Y+ e: L- B6 e3 T) m9 ?
    for j=1:num_inter
: M% }! W- Y- D2 t0 u9 C       w2=Py_Pgama(d,wpeak,wframe,1,sr);   % 先进行Py投影和 Pgama投影7 Q1 F+ V3 A5 ~. Y# `' q7 A1 F
       w0=iswt(pswa,w2,Lo_R,Hi_R);         % 再进行Pv投影9 Q: p$ }) S- V: g2 Q6 f* d
       [a,d]=swt(w0,level,Lo_D,Hi_D);      % Pv8 G* a5 X* _, e1 i) Y5 X
    end5 s3 _' b7 f/ O$ h7 z
     pswa=iswt(swa(level,:),w2,Lo_R,Hi_R); % 计算重建信号   
/ R( V* g: ~! u: s/ k     " b! ?$ r/ `" s( d/ A$ v0 L6 b
xcrr=y'-pswa; % 重建误差
: D  k1 [& e" D. i7 ^. J5 C4 xfigure;%figure3
( X! z6 o5 V- L2 E! Wsubplot(411)
8 G. P) r! T( q- T4 o0 V2 I1 }( @plot(y(1:points),'r');
* j6 |7 F; z  h! uaxis tight;
2 p2 R8 g1 D2 ?# L9 L( dsubplot(412)
, p' w! @/ X1 X" t; I7 {plot(signal(1:points),'r');
/ n- r6 L) Q2 g& k: qaxis tight;
; K" e! x: g/ r+ Usubplot(413)0 Y; j' p: ~0 x+ r1 F/ B( {* Q
plot(pswa(1:points)); 1 o9 `, D7 `1 u: d+ F, d5 q: C
axis tight;3 y/ N; N7 ^/ i% r! j
subplot(414)0 F: B: H! A9 ^& J8 q  }- ]
plot(xcrr(1:points));
- G+ F7 O( w  l1 Waxis tight;  a/ {' K. w. c& D
figure%figure4) C1 x) J5 B9 u# _5 r" F5 ^
subplot(611)
$ n% r7 o, H- }: s# r% m! l7 ~plot(signal);
, c1 v# @  U+ {9 usubplot(612)
, _& a: D1 W$ ?/ I/ ~plot(D5_wpeak);2 [" H: W, O# X' z2 J1 ^8 q7 g  h
subplot(613)5 G9 S8 u3 E6 {( i+ q5 u
plot(D4_wpeak);
% M: b# ]6 c6 ^8 M% S! s) N1 ?1 nsubplot(614)/ ~! E& ^% h3 C8 g" h
plot(D3_wpeak);
" M8 S8 F/ n' `6 K; y6 ]/ }subplot(615)' M5 s- L  n* o# t% g' g1 D
plot(D2_wpeak);
1 d% N9 x, h$ Q3 [1 W( w2 esubplot(616)3 {. Y7 }8 e1 g  _( |
plot(D1_wpeak);
/ |8 C" s% ~0 p8 O[value,index]=max(abs(D2_wpeak))9 D: L0 X5 s  m; y# m+ E) V
[value,index]=max(abs(D1_wpeak))
更多图片 小图 大图
组图打开中,请稍候......
"真诚赞赏,手留余香"
还没有人打赏,支持一下
楼主热帖
帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
  • TA的每日心情
    郁闷
    2017-6-11 22:50
  • 签到天数: 4 天

    连续签到: 2 天

    [LV.2]偶尔看看I

    累计签到:5 天
    连续签到:1 天
    发表于 2009-5-12 14:43:18 | 显示全部楼层
    先用MATLAB自带的小波工具箱选择不同小波函数试试效果,你这个分解没有产生奇异点啊
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】

    该用户从未签到

    尚未签到

    发表于 2010-3-31 20:39:44 | 显示全部楼层
    可以找到突变点吗?我想了解了解,能不能说一下啊
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】

    该用户从未签到

    尚未签到

    发表于 2010-4-21 10:15:30 | 显示全部楼层
    回复 3# 阳光zzl ( \7 y- v' b$ u' T4 m4 l

    9 y* v7 h' ~3 S, Q/ W. ~! D
    / G' v0 l5 T( |& y2 @# m  I    我也来请教~·
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】

    该用户从未签到

    尚未签到

    发表于 2010-7-21 11:47:20 | 显示全部楼层
    我也对模极大值感兴趣
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】

    该用户从未签到

    尚未签到

    发表于 2010-8-8 19:00:47 | 显示全部楼层
    对小波的模极大值,也需要学习的~有些不明白~
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
    您需要登录后才可以回帖 登录 | 立即加入

    本版积分规则

    招聘斑竹

    小黑屋|手机版|APP下载(beta)|Archiver|电力研学网 ( 赣ICP备12000811号-1|赣公网安备36040302000210号 )|网站地图

    GMT+8, 2026-3-19 08:13

    Powered by Discuz! X3.5 Licensed

    © 2001-2025 Discuz! Team.

    快速回复 返回顶部 返回列表