给大家分享我自己编的程序-连续小波变换

pengzk · 发表于 2007-4-19 04:20

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

Function [WT, FreqBins, Scales] = 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;

pengzk · 发表于 2007-4-19 19:58

Function [WT, FreqBins, Scales] = CWT_Morlet(Sig, WinLen, nLevel);

改为

Function [WT, FreqBins, Scales] = Wavelet_Morlet(Sig, WinLen, nLevel);

eight · 发表于 2007-4-19 20:23

原帖由 pengzk 于 2007-4-19 19:58 发表
Function = CWT_Morlet(Sig, WinLen, nLevel);

改为

Function = Wavelet_Morlet(Sig, WinLen, nLevel);

你编辑一下自己的帖子（进入帖子后，在1楼的右下角点击“编辑”便可），这样方便大家下载

wuhanxy123 · 发表于 2007-7-31 21:13

该程序有问题: 首先从if到第一个end之间的语句与程序无关, 其次没有实现时间方向的平移动, 每个频率水平只求了一个时窗的小波系数!!

zhlong · 发表于 2007-7-31 21:24

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是一个向量，已经是实现了时间方向的平移。

wxr1000 · 发表于 2007-7-31 21:26

为什么要自己编程序呢？matlab中不是有现成的吗？

wuhanxy123 · 发表于 2007-7-31 22:39

晕, 如果m固定, 则频率参数a固定, 时间向量T的长固定, 这个长度不等于整体信号的长度, 所以必须还要移动窗口,使之覆盖整个信号, 非再加循环不可.

[ 本帖最后由 zhlong 于 2007-8-1 07:16 编辑 ]

zhlong · 发表于 2007-8-1 08:08

楼主还发了这个函数的应用实例，你可以验证一下。
http://forum.vibunion.com/forum/ ... highlight=%2Bpengzk

zhlong · 发表于 2007-8-1 10:44

楼主的程序个人觉得没有问题，你可以参考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;
[WT, FreqBins, Scales] = CWT_Morlet(sig,WinLen,512);
FreqBins = FreqBins * SampFreq;
clf
set(gcf,'Position',[20 100 300 220]);
set(gcf,'Color','w');
pcolor(t,FreqBins,abs(WT));
colormap jet;
shading interp;
axis([min(t) max(t) min(FreqBins) max(FreqBins)]);
colorbar;
ylabel('Frequency / Hz');
xlabel('Time / sec');

[ 本帖最后由 zhlong 于 2007-8-1 10:45 编辑 ]

zhlong · 发表于 2007-11-12 13:51

FreqBins = logspace(log10(fmin),log10(0.5),nLevel); % 将频率轴在分析范围内等对数坐标划分
Scales = fmax*ones(size(FreqBins))./ FreqBins; % 计算响应的尺度参数

不知道这个频率和尺度之间的关系是根据什么原理确定的？？

破凰 · 发表于 2007-11-18 16:49

原帖由 zhlong 于 2007-11-12 13:51 发表

 登录/注册后可看大图

不知道这个频率和尺度之间的关系是根据什么原理确定的？？

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 [WT, FreqBins, Scales] = 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:11

很好!

[ 本帖最后由 kevin19821 于 2007-12-18 11:13 编辑 ]

bitww · 发表于 2007-12-20 21:52

受益非浅呀

hyl2323 · 发表于 2007-12-21 08:02

好东西，我回去研究一下。

kevin19821 · 发表于 2007-12-21 15:48

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是做啥的?

账号		自动登录	找回密码
密码			我要加入

[小波] 给大家分享我自己编的程序-连续小波变换

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

评分

更正

回复 #7 wuhanxy123 的帖子

评分