在文章“給大家分享我自己編的程式-連續小波變換” 中，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.