[转载]论坛中某位达人自己编写的Morlet连续小波变换程序

原文地址：论坛中某位达人自己编写的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小波变换的函数为：

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的值即可

以上程序转自振动联盟论坛：http://www.chinavib.com/thread-40900-1-2.html

实际上，matlab有自己的Morlet小波函数，即在小波分析的工具箱中自带着，帖这个程序的目的，是为了弄清楚小波变换的原理