在帖子“给大家分享我自己编的程序-连续小波变换” 中，pengzk版友给出了morlet小波变换的源代码，但其中的许多参数和语句意义不够明确，这就给一些希望了解连续小波变换实现方法的版友带来不便。因此，本帖将对连续小波变换的实现原理做个小结，希望对各位有所帮助。
首先说明的是，在Matlab的小波工具箱和pengzk版友提供的程序中，连续小波变换都是依据以下原理实现的：连续小波变换可以看成是原信号与小波基进行卷积的结果。
下面我以自编的连续morlet小波变换程序为例说明利用卷积方法实现连续小波变换的过程(参见程序注释)。其中，所用morlet小波的定义为
程序如下：
function wcoefs = mymorletcwt(Sig,Scales,fc,fb)
%==================%
% Continuous Wavelet Transform using Morlet function
%======输入======%
% Sig: 输入信号
% Scales: 输入的尺度序列
% fc: morlet小波中心频率 (默认为2)
% fb: morlet小波带宽参数 (默认为2)
%======输出======%
% wcoefs: morlet小波变换计算结果
%==================%
if (nargin <= 1),
error('At least 2 parameters required');
end;
if (nargin < 4),
fb = 2;
elseif (nargin < 3),
fc = 2;
end;
wavsupport=8;
% 默认morlet小波的支撑区为[-8,8]
nLevel=length(Scales);
% 尺度的数目
SigLen = length(Sig);
% 信号的长度
wcoefs = zeros(nLevel, SigLen);
% 分配计算结果的存储单元
for m = 1:nLevel
% 计算各尺度上的小波系数
a = Scales(m);
% 提取尺度参数
t = -round(a*wavsupport):1:round(a*wavsupport);
% 在尺度a的伸缩作用下，此时小波函数的支撑区会变为[-a*wavsup,a*wavsup]，采样频率为1Hz
Morl = fliplr((pi*fb)^(-1/2)*exp(i*2*pi*fc*t/a).*exp(-t.^2/(fb*a^2)));
% 计算当前尺度下的小波函数，按小波变换的定义这里需要倒置
temp = conv(Sig,Morl) / sqrt(a);
% 计算信号与当前尺度下小波函数的卷积
d=(length(temp)-SigLen)/2;
% 由于卷积计算所得结果的长度可能远远大于原信号，只需提取按原信号的长度提取中间部分的系数
first = 1+floor(d);
% 区间的起点
wcoefs(m,:)=temp(first:first+SigLen-1);
end
从以上程序可以看出，基于卷积原理的连续小波变换实现的关键是求出某一尺度下的小波函数离散化后的序列。该序列可以通过对该尺度下的小波函数进行采样求得，采样区间为此时小波函数的支撑区，采样频率可取为1Hz。
注：(1)按我的理解，采样频率取得越大，计算结果也越精确，但我在测试中发现，采样频率取得太高反而会影响分析结果的精度，在本例中采样频率最好取为1Hz。
(2)在小波工具箱的cwt函数中，某一尺度下的小波函数采样序列是直接对母小波采样序列伸缩而得。
下面利用zhlong给出的例子对mymorletcwt进行测试，并与小波工具箱中cwt所得结果进行对比。
clc;
clear;
SampFreq = 30;
t=1/SampFreq: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));
fmax = 0.5;
% 最高分析频率(归一化频率)
fmin = 0.005;
% 最低分析频率(归一化频率)
fb = 4 ;
% 取cmor4-2小波进行实验，带宽参数为4
fc = 2;
% 中心频率2Hz
totalscal = 512;
% 所取尺度的数目
FreqBins = linspace(fmin,fmax,totalscal);
% 将频率轴在分析范围内等间隔划分
Scales = fc./ FreqBins;
% 计算相应的尺度参数
wcoefs = mymorletcwt(sig,Scales,fc,fb);
RealFreqBins = FreqBins * SampFreq;
% 尺度所对应的实际频率
%====本帖方法的结果====%
figure(1)
pcolor(t,RealFreqBins,abs(wcoefs));
colormap jet;
shading interp;
axis([min(t) max(t) min(RealFreqBins) max(RealFreqBins)]);
colorbar;
ylabel('Frequency / Hz');
xlabel('Time / sec');
%==小波工具箱中的cwt结果==%
figure(2)
MWT=cwt(sig,Scales,'cmor4-2');
pcolor(t,RealFreqBins,abs(MWT));
colormap jet;
shading interp;
colorbar;
axis([min(t) max(t) min(RealFreqBins) max(RealFreqBins)]);
ylabel('Frequency / Hz');
xlabel('Time / sec');
测试结果如下：
本帖方法的结果
小波工具箱中的cwt结果
对比两图可见，两种方法的精度相仿。我通过查看cwt的源代码发现，两者的区别在于小波工具箱中的cwt是按先对小波函数积分，卷积后再微分的方法来求小波系数的。至于这样做的根据就不得而知了。对于其它小波，也应该可以按以上原理实现，特别是那些可以用解析形式表示的小波，如Gaussian小波、Shannon小波、B样条小波等。
小结：按卷积原理实现的连续小波变换的方法简单直观，并且个人认为其本质上是一种矩形数值积分法。但这种方法的计算精度和速度可能不如其它方法，如更高精度的数值积分法、调频Z变换法、梅林变换法等。在此，偶也希望其它版友能积极发帖对这几种方法进行讨论。
如何利用Matlab进行小波分析
小波分析(wavelet analysis), 或小波转换(wavelet transform)是指用有限长或快速衰减的、称为母小波(mother wavelet)的振荡波形来表示信号。该波形被缩放和平移以匹配输入的信号。
小波变换的概念是由法国从事石油信号处理的工程师J.Morlet在1974年首先提出的，通过物理的直观和信号处理的实际经验的需要建立了反演公式，当时未能得到数学家的认可。正如1807年法国的热学工程师J.B.J.Fourier提出任一函数都能展开成三角函数的无穷级数的创新概念未能得到著名数学家J.L.Lagrange，P.S.Laplace以及A.M.Legendre的认可一样。幸运的是，早在七十年代，A.Calderon表示定理的发现、Hardy空间的原子分解和无条件基的深入研究为小波变换的诞生做了理论上的准备，而且J.O.Stromberg还构造了历史上非常类似于当前的小波基;1986年著名数学家Y.Meyer偶然构造出一个真正的小波基，并与S.Mallat合作建立了构造小波基的统一方法加多尺度分析之后，小波分析才开始蓬勃发展起来，其中比利时女数学家I.Daubechies撰写的《小波十讲(Ten Lectures on Wavelets)》对小波的普及起了重要的推动作用。它与Fourier变换、窗口Fourier变换(Gabor变换)相比，这是一个时间和频率的局域变换，因而能有效的从信号中提取信息，通过伸缩和平移等运算功能对函数或信号进行多尺度细化分析(MultiscaleAnalysis)，解决了Fourier变换不能解决的许多困难问题，从而小波变化被誉为"数学显微镜"，它是调和分析发展史上里程碑式的进展。
本文以《Wavelet transforms and their applications toturbulence》和《A practicalguide to wavelet analysis》两篇文章的研究为基础，利用MATLAB，结合西北地区ET0,研究分析了其周期变化。并给出具体代码：
%WAVETEST Example Matlab script for WAVELET, using NINO3 SST dataset
%
% See "http://paos.colorado.edu/research/wavelets/"
% Written January 1998 by C. Torrence
%
% Modified Oct 1999, changed Global Wavelet Spectrum (GWS) to be sideways,
% changed all "log" to "log2", changed logarithmic axis on GWS to
% a normal axis.
load 'sst.txt' % input SST time series
sst = sst;
%------------------------------------------------------ Computation
% normalize by standard deviation (not necessary, but makes it easier
% to compare with plot on Interactive Wavelet page, at
% "http://paos.colorado.edu/research/wavelets/plot/"
variance = std(sst)^2;
sst = (sst - mean(sst))/sqrt(variance) ;
n = length(sst);
dt = 1 ;
time = [0:length(sst)-1]*dt + 1956.0 ; % construct time array
xlim = [1956,2011]; % plotting range
pad = 1; % pad the time series with zeroes (recommended)
dj = 0.25; % this will do 4 sub-octaves per octave
s0 = 2*dt; % this says start at a scale of 6 months
j1 = 7/dj; % this says do 7 powers-of-two with dj sub-octaves each
lag1 = 0.72; % lag-1 autocorrelation for red noise background
mother = 'Morlet';
% Wavelet transform:
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);
power = (abs(wave)).^2 ; % compute wavelet power spectrum
% Significance levels: (variance=1 for the normalized SST)
[signif,fft_theor] = wave_signif(1.0,dt,scale,0,lag1,-1,-1,mother);
sig95 = (signif')*(ones(1,n)); % expand signif --> (J+1)x(N) array
sig95 = power ./ sig95; % where ratio > 1, power is significant
% Global wavelet spectrum & significance levels:
global_ws = variance*(sum(power')/n); % time-average over all times
dof = n - scale; % the -scale corrects for padding at edges
global_signif = wave_signif(variance,dt,scale,1,lag1,-1,dof,mother);
% Scale-average between El Nino periods of 2--8 years
avg = find((scale >= 2) & (scale < 8));
Cdelta = 0.776; % this is for the MORLET wavelet
scale_avg = (scale')*(ones(1,n)); % expand scale --> (J+1)x(N) array
scale_avg = power ./ scale_avg; % [Eqn(24)]
scale_avg = variance*dj*dt/Cdelta*sum(scale_avg(avg,:)); % [Eqn(24)]
scaleavg_signif = wave_signif(variance,dt,scale,2,lag1,-1,[2,7.9],mother);
whos
%------------------------------------------------------ Plotting
%--- Contour plot wavelet power spectrum
subplot('position',[0.1 0.37 0.65 0.28])
levels = [0.0625,0.125,0.25,0.5,1,2,4,8,16] ;
Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
contour(time,log2(period),log2(power),log2(levels)); %*** or use 'contourfill'
%imagesc(time,log2(period),log2(power)); %*** uncomment for 'image' plot
xlabel('Time (year)')
ylabel('Period (years)')
title('Wavelet Power Spectrum')
set(gca,'XLim',xlim(:))
set(gca,'YLim',log2([min(period),max(period)]), ...
'YDir','reverse', ...
'YTick',log2(Yticks(:)), ...
'YTickLabel',Yticks)
% 95% significance contour, levels at -99 (fake) and 1 (95% signif)
hold on
contour(time,log2(period),sig95,[-99,1],'k');
hold on
% cone-of-influence, anything "below" is dubious
plot(time,log2(coi),'k')
hold off
具体分析：西北地区的年均ET0存在2-3a的显著性震荡周期和6a的准周期震荡。对于2-3a的显著周期震荡，分别在20世纪60年代初期到60年代末期和20世纪80年代末期到21世纪初期较为显著，并且处于低值期，2005年以后进入减小期；对于6a的准周期震荡而言，分别在20世纪60年代后期到20世纪70年代前期、20世纪70年代后期到20世纪80年代初期和20世纪80年代中期到21世纪初期，并且目前处于自2006年以来ET0的低值区。值得注意的是，在超过28年的分析中，等值线几乎全为红色，且没有闭合的中心，这主要是由于ET0的时间序列仅为56年，超过28年的周期不能明显地表示出来。
参考文献
[1] Marie Farge. Wavelet transforms and their applications to turbulence[J]. Annual Review of Fluid Mechanics, 1992,24: 395-457.
[2] Christopher Torrence, Gilbert P. Compo. A practical guide to wavelet analysis[J]. Bulletin of the American Meteorological Society, 1998,79: 61-78.