给大家分享我自己编的程序-连续小波变换
Function = CWT_Morlet(Sig, WinLen, nLevel);%============================================================%
%Continuous Wavelet Transform using Morlet function
% Sig : 信号
% WinLen : 小波函数在尺度参数a=1时的长度 （默认为 10）
% nLevel : 频率轴划分区间数 （默认为1024）
%
% WT:返回的小波变换计算结果
%FreqBins :返回频率轴划分结果（归一化频率，最高频率为0.5）
% Scales: 返回与频率轴划分值相对应的尺度划分 （频率0.5对应的尺度为1）
%============================================================%
if (nargin == 0),
error('At least 1 parameter required');
end;
if (nargin < 4),
iShow = 1;
elseif (nargin < 3),
nLevel = 1024;
elseif (nargin < 2),
WinLen = 10;
end;
Sig = hilbert(real(Sig)); % 计算信号的解析信号
SigLen = length(Sig); % 获取信号的长度
fmax = 0.5; % 设置最高分析频率
fmin = 0.005; % 设置最低分析频率
FreqBins = logspace(log10(fmin),log10(0.5),nLevel); % 将频率轴在分析范围内等对数坐标划分
Scales = fmax*ones(size(FreqBins))./ FreqBins; % 计算响应的尺度参数
omg0 = WinLen / 6; % 按给定的小波长度计算相应的小波函数中心频率
WT = zeros(nLevel, SigLen); % 分配计算结果的存储单元
wait = waitbar(0,'Under calculation, please wait...');
for m = 1:nLevel,
waitbar(m/nLevel,wait);
a = Scales(m); % 提取尺度参数
t = -round(a*WinLen):1:round(a*WinLen);
Morl = pi^(-1/4)*exp(i*2*pi*0.5*t/a).*exp(-t.^2/2/(2*omg0*a)^2); % 计算当前尺度下的小波函数
temp = conv(Sig,Morl) / sqrt(a); % 计算信号与小波函数的卷积
WT(m,:) = temp(round(a*WinLen)+1:length(temp)-round(a*WinLen));
end;
close(wait);
WT = WT / WinLen;
更正
Function = CWT_Morlet(Sig, WinLen, nLevel);改为
Function = Wavelet_Morlet(Sig, WinLen, nLevel); 原帖由 pengzk 于 2007-4-19 19:58 发表
Function= CWT_Morlet(Sig, WinLen, nLevel);
改为
Function= Wavelet_Morlet(Sig, WinLen, nLevel);
你编辑一下自己的帖子（进入帖子后，在1楼的右下角点击“编辑”便可），这样方便大家下载 该程序有问题: 首先从if到第一个end之间的语句与程序无关, 其次没有实现时间方向的平移动, 每个频率水平只求了一个时窗的小波系数!! if (nargin == 0),
error('At least 1 parameter required');
end;
这个用于检查函数的输入。楼上说的是这个判断语句吗？
t = -round(a*WinLen):1:round(a*WinLen);
Morl = pi^(-1/4)*exp(i*2*pi*0.5*t/a).*exp(-t.^2/2/(2*omg0*a)^2); % 计算当前尺度下的小波函数
这里的t是一个向量，已经是实现了时间方向的平移。 为什么要自己编程序呢？matlab中不是有现成的吗？ 晕, 如果m固定, 则频率参数a固定, 时间向量T的长固定, 这个长度不等于整体信号的长度, 所以必须还要移动窗口,使之覆盖整个信号, 非再加循环不可.
[ 本帖最后由 zhlong 于 2007-8-1 07:16 编辑 ]
回复 #7 wuhanxy123 的帖子
楼主还发了 这个函数的应用实例，你可以验证一下。http://forum.vibunion.com/forum/viewthread.php?tid=41381&highlight=%2Bpengzk 楼主的程序个人觉得没有问题，你可以参考tftb时频工具箱中的tfrscalo函数，比较一下这两个程序。
SampFreq = 30;
t=0:1/SampFreq:4;
sig = sin(12*pi*t);
sig(1:end/2) = sig(1:end/2) + sin(6*pi*t(1:end/2));
sig(end/2+1:end) = sig(end/2+1:end) + sin(18*pi*t(end/2+1:end));
WinLen = 10;
= CWT_Morlet(sig,WinLen,512);
FreqBins = FreqBins * SampFreq;
clf
set(gcf,'Position',);
set(gcf,'Color','w');
pcolor(t,FreqBins,abs(WT));
colormap jet;
shading interp;
axis();
colorbar;
ylabel('Frequency / Hz');
xlabel('Time / sec');
[ 本帖最后由 zhlong 于 2007-8-1 10:45 编辑 ] FreqBins = logspace(log10(fmin),log10(0.5),nLevel); % 将频率轴在分析范围内等对数坐标划分
Scales = fmax*ones(size(FreqBins))./ FreqBins; % 计算响应的尺度参数
不知道这个频率和尺度之间的关系是根据什么原理确定的？？ 原帖由 zhlong 于 2007-11-12 13:51 发表 http://www.chinavib.com/forum/images/common/back.gif
不知道这个频率和尺度之间的关系是根据什么原理确定的？？
lz好像把morlet小波的定义整错了。
我看的书上的morlet小波的定义是
其中，fc是小波中心频率，fb是带宽参数。在lz贴的程序中，中心频率为omg0 = WinLen / 6;而在计算时小波函数用的是
Morl = pi^(-1/4)*exp(i*2*pi*0.5*t/a).*exp(-t.^2/2/(2*omg0*a)^2) （1）
与定义相比可以看出，程序里注释的中心频率和带宽参数弄反了，所用morlet小波的中心频率应是0.5Hz。然后按照尺度和频率关系（小波中心频率除以尺度等于归一化频率）就有
Scales = 0.5*ones(size(FreqBins))./ FreqBins
程序里fmax刚好又等于0.5，所以图上的标注是对的，但这只是种偶然，如果将式（1）中的0.5改为其它数，得到的频率坐标就不对。
PS：个人认为lz提供的函数的输入参数WinLen不仅与母小波的长度（支撑区）有关，而且与morlet小波的带宽参数有关（通过WinLen得到了omg0）。
我把这个函数修改了一下
function = CWT_Morlet(Sig, WinLen, nLevel);
%============================================================%
%Continuous Wavelet Transform using Morlet function
% Sig : 信号
% WinLen : 母小波离散化后的长度的一半 （默认为 10）
% nLevel : 频率轴划分区间数 （默认为1024）
%
% WT:返回的小波变换计算结果
%FreqBins :返回频率轴划分结果（归一化频率，最高频率为0.5）
% Scales: 返回与频率轴划分值相对应的尺度划分 （频率0.5对应的尺度为1）
%============================================================%
if (nargin == 0),
error('At least 1 parameter required');
end;
if (nargin < 4),
iShow = 1;
elseif (nargin < 3),
nLevel = 1024;
elseif (nargin < 2),
WinLen = 10;
end;
Sig = hilbert(real(Sig)); % 计算信号的解析信号
SigLen = length(Sig); % 获取信号的长度
fmax = 0.5; % 设置最高分析频率
fmin = 0.005; % 设置最低分析频率
fc = 0.5; %小波中心频率
omg0 = WinLen / 6; % 小波带宽参数
FreqBins = logspace(log10(fmin),log10(0.5),nLevel); % 将频率轴在分析范围内等对数坐标划分
Scales = fc*ones(size(FreqBins))./ FreqBins; % 计算相应的尺度参数
WT = zeros(nLevel, SigLen); % 分配计算结果的存储单元
wait = waitbar(0,'Under calculation, please wait...');
for m = 1:nLevel,
waitbar(m/nLevel,wait);
a = Scales(m); % 提取尺度参数
t = -round(a*WinLen):1:round(a*WinLen);
Morl = pi^(-1/4)*exp(i*2*pi*fc*t/a).*exp(-t.^2/2/(2*omg0*a)^2); % 计算当前尺度下的小波函数
temp = conv(Sig,Morl) / sqrt(a); % 计算信号与小波函数的卷积
WT(m,:) = temp(round(a*WinLen)+1:length(temp)-round(a*WinLen));
end;
close(wait);
WT = WT / WinLen;
如果要改变中心频率，直接改变fc的值即可。
[ 本帖最后由 zhlong 于 2007-12-21 20:06 编辑 ] 很好!
[ 本帖最后由 kevin19821 于 2007-12-18 11:13 编辑 ] 受益非浅呀 好东西，我回去研究一下。 Morl = pi^(-1/4)*exp(i*2*pi*0.5*t/a).*exp(-t.^2/2/(2*omg0*a)^2);
exp(-t.^2/2/(2*omg0*a)^2); 中的(2*omg0*a)^2是做啥的?