|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 guiru8889 于 2011-5-26 15:37 编辑
RWT的调用格式是:
RWT(x,nvoice,wavelet,oct,scale)
x:待分析信号
nvoice:两个尺度之间的取nvoice个间隔点
wavelet:小波名字
oct:原文中对oct的引用是noctave = floor(log2(n))-oct;看来oct是控制尺度范围的,但是尺度范围也由信号长度控制即: floor(log2(n))
scale:起始尺度值,后面将以2倍增,每两倍之间插入nvoice个间隔点。体现代码:
for jo = 1:noctave,
for jv = 1:nvoice,
qscale = scale .* (2^(jv/nvoice));
end
scale = scale .*2;
end
RWT的原理:
RWT采用的是循环卷积的形式来计算连续小波变换,当然只能是近似了。假设信号的长度为n,RWT也对相关小波函数采样n个点,当然采样的步长就和尺度有关系了。
本来采用循环卷积就应当直接求小波函数的采样值,然后求其fft,再和信号的fft相乘,最后求ifft,但是RWT为了节省一步计算量,直接对小波函数的傅立叶变换采样,具体体现是:
xi = [ (0: (n/2)) (((-n/2)+1):-1) ] .* (2*pi/n);
omega = n .* xi ./ qscale ;
if strcmp(wavelet,'Gauss'),
window = exp(-omega.^2 ./2);
elseif strcmp(wavelet,'DerGauss'),
window = i.*omega.*exp(-omega.^2 ./2);
elseif strcmp(wavelet,'Sombrero'),
window = (omega.^2) .* exp(-omega.^2 ./2);
elseif strcmp(wavelet,'Morlet'),
window = exp(-(omega - omega0).^2 ./2) - exp(-(omega.^2 + omega0.^2)/2);
end
xi我认为就是在一个2pi周期内取等间隔的n个值,omega = n .* xi ./ qscale ;把频率和尺度联系了起来(备注:omega = n .* xi ./ qscale ,中的n其实实不应该有,因为这里的n只是分析数据的长度,所以它影响的是尺度值,也就是qscale,如果不乘以这个n,在调用rwt函数的时候,尺度范围的输入值就应当减小n倍,你可以修改一下试试。)
接下来的事情就是ifft了。不多说了。
不过出于好奇心,我实现了一个对小波函数直接采样的程序,具体就是:
xi = [ (0: (n/2)) (((-n/2)+1):-1) ];
omega = xi ./qscale ;
if strcmp(wavelet,'Gauss'),
window = 1/sqrt(2*pi).*exp(-omega.^2 ./2);
elseif strcmp(wavelet,'DerGauss'),
window = 1/sqrt(2*pi).*(-omega).*exp(-omega.^2 ./2);
end
我只写了高斯函数和其一阶导数的采样,其他的没写。程序运行结果和源程序一样。注意尺度的取值范围很关键。
源程序的运行结果
源程序的运行结果
我修改之后的程序的运行结果
|
评分
-
1
查看全部评分
-
|