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

總結

以上是生活随笔為你收集整理的matlab 小波变换_连续小波变换实现方法的总结及其程序详解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。