WAVELAB中的RWT连续小波变换分析
本帖最后由 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
我只写了高斯函数和其一阶导数的采样，其他的没写。程序运行结果和源程序一样。注意尺度的取值范围很关键。
源程序的运行结果
我修改之后的程序的运行结果
非常感谢楼主的解读啊 本帖最后由 甜心宝贝 于 2011-6-7 19:31 编辑
楼主我比较笨，看了许久终于是弄懂了，不过还是有些疑问，对于楼主所说的omega =n .* xi ./ qscale ，中的n其实实不应该有，这句没明白，程序中
xi = [ (0: (n/2)) (((-n/2)+1):-1) ] .* (2*pi/n);
omega =n .* xi ./ qscale ;这不是和楼主下面的是一个意思么，除以了n下面还是乘以了n，如果omega =n .* xi ./ qscale 中不乘以n的话那尺度就会增大n，所以这个程序中还是应该乘以n的吧，我感觉楼主想要表达的应该也是这个意思，还有楼主你实现了的一个对小波函数直接采样的程序我还是有些没懂，为什么window = 1/sqrt(2*pi).*exp(-omega.^2 ./2)这句中要除以1/sqrt(2*pi）呢？望不吝赐教，我比较笨，学习中，先谢谢了
不乘以n，输入的是真实尺度范围，乘以n之后，需要输入的尺度就不是真实的了，是扩大n倍的尺幅范围了。
window = 1/sqrt(2*pi).*exp(-omega.^2 ./2)这句中之所以要除以1/sqrt(2*pi）是因为高斯小波函数本来就需要的。
页:
[1]