有做小波脊线提取的嘛？

avaxu

最近在做小波脊线提取，没有一点思路，有相应的文献或者程序可以借鉴不？

看到有从模最大值和相位方面2种提取方

破凰

原帖由 avaxu 于 2007-6-15 11:06 发表

                               
登录/注册后可看大图

最近在做小波脊线提取，没有一点思路，有相应的文献或者程序可以借鉴不？

看到有从模最大值和相位方面2种提取方

模最大值法好像不是求小波脊线的吧，它是另一种求瞬时频率的方法。小波脊线提取一般都是采用一种迭代算法。介绍该算法的文献很多，例如“基于小波脊线的刮水器电动机机械特性测试”。

avaxu

谢谢，破凰！  一直将小波脊线理解为最大值的连线，就认为是模最大值。

gemboy

Multiridge Detection and Time–Frequency Reconstruction.pdf

zhlong

一篇硕士学位论文看到如下叙述，不知道是不是叙述有误？

batizhang

好像不止一篇文章是这么说的，哪位大侠出来说明一下，这种说法对不对啊

破凰

小波脊线上的点对应的是该时刻小波系数模极大值的点，是我弄错了。

yao1983

clc
clear
%设置a的初值；
a(1)=0.02;
Te=0.001;
t=0:0.001:2;
s=cos(2*pi*30*(t.^2));
s=hilbert(s);
hudu=4*pi;
%用morlet小波进行分析；
kk=length(t);
k=0;
aa=0;
bb=0;
for j=1:1000
wa=cwt(s,a(1),'morl');
fai=angle(wa);
for i=2:kk
a(i)=hudu/((fai(i)-fai(i-1))/Te); %求Db；
if abs(a(i)-a(i-1))/(abs(a(i-1)))<0.001 %判断是否满足条件；
k=k+1;
aa=a(i);
break;
end
end
ar(k)=abs(aa);%赋值；
a(1)=ar(k);%继续求a;
end
plot(ar);

不知道哪点出现错误了,请高手指点

破凰

（1）用复的morlet小波
（2）求出相位后要解卷绕
（3）求小波脊线的公式不对，再仔细看下资料。

yao1983

楼上的高手:
你的第二个意见:（2）求出相位后要解卷绕 是 什么意思呀,还有（3）求小波脊线的公式不对,我是按照你提供的文章写的程序,公式不可能不对呀,还望你再次指点一下

破凰

解卷绕用unwrap。
a(i)=hudu/((fai(i)-fai(i-1))/Te); 这个公式中多除了个采样频率。当然如果你最后要转化为瞬时频率的话，也可以在这里先除，但此时的小波脊线不是反映的真实尺度。

yao1983

高手我按照你说的改了一下程序，可是结果还是不对，还的麻烦你一下，再帮我看看程序，我是新手，对小波也是刚刚入门
clc
clear
%设置a的初值；
a(1)=0.02;
Te=0.001;
t=0:0.001:2;
s=cos(2*pi*30*(t.^2));
s=hilbert(s);
%hudu=4*pi;
%用morlet小波进行分析；
kk=length(t);
k=0;
aa=0;
bb=0;
for j=1:1000
wa=cwt(s,a(1),'cmor2-4*pi');
fai=angle(wa);
fai=unwrap(fai);
for i=2:kk
a(i)=4*pi/((fai(i)-fai(i-1))/Te); %求Db；
if abs(a(i)-a(i-1))/(abs(a(i-1)))<0.001 %判断是否满足条件；
k=k+1;
aa=a(i);
break;
end
end
ar(k)=abs(aa);%赋值；
a(1)=ar(k);%继续求a;
end
plot(ar);
下面是程序的输出结果

破凰

刚刚又看了一遍，发现你把算法理解错了。

我编的程序为：

clc;
clear all;
a(1)=25;
Te=0.001;
t=0:0.001:1;
s=cos(2*pi*(30*(t.^2)+40*t));  %改了一下原信号
%s=hilbert(s);
%hudu=4*pi;
%用morlet小波进行分析；
kk=length(t);
k=0;
aa=0;
bb=0;
w0=2*pi;  %cmor2-1小波的中心频率为2pi
for j=2:kk
    while 1
        wa=cwt(s,a(1),'cmor2-1');   
        fai=angle(wa);
        fai=unwrap(fai);
        dfa=abs(fai(j)-fai(j-1));
        a0=w0*1/dfa;
        if abs((a0-a(1))/a(1))<0.02  %阈值不宜过小
            ar(j-1)=a(1);
            a(1)=a0;
            break;
        end
        a(1)=a0;
    end
end
f=1./(ar*Te);
plot(t(1:kk-1),f);

结果

没怎么写注释，你对着算法再好好看下吗。

[ 本帖最后由 破凰 于 2007-11-23 20:53 编辑 ]

yao1983

多谢多谢！:handshake

liuhui1983

在凤凰编的程序的基础上稍加改进,不过仍旧存在很大的边界问题,还请高手指点指点!!
clc;
clear all;
a(1)=25;
Te=0.001;
t=0:0.001:1;
s=cos(2*pi*(30*(t.^2)+40*t));  %改了一下原信号
s=hilbert(s);
%hudu=4*pi;
%用morlet小波进行分析；
kk=length(t);
k=0;
aa=0;
bb=0;
w0=2*pi;  %cmor2-1小波的中心频率为2pi
for j=2:kk
    while 1
        wa=cwt(s,a(1),'morl');   
        fai=angle(wa);
        fai=unwrap(fai);
        dfa=abs(fai(j)-fai(j-1));
        a0=w0*1/dfa;
        if abs((a0-a(1))/a(1))<0.9  %阈值不宜过小
            ar(j-1)=a(1);
            a(1)=a0;
            break;
        end
        a(1)=a0;
    end
end
f=1./(ar*Te);
plot(t(1:kk-1),f);