在帖子“給大家分享我自己編的程序-連續(xù)小波變換” 中，pengzk版友給出了morlet小波變換的源代碼，但其中的許多參數(shù)和語句意義不夠明確，這就給一些希望了解連續(xù)小波變換實現(xiàn)方法的版友帶來不便。因此，本帖將對連續(xù)小波變換的實現(xiàn)原理做個小結(jié)，希望對各位有所幫助。

首先說明的是，在Matlab的小波工具箱和pengzk版友提供的程序中，連續(xù)小波變換都是依據(jù)以下原理實現(xiàn)的：連續(xù)小波變換可以看成是原信號與小波基進(jìn)行卷積的結(jié)果。

下面我以自編的連續(xù)morlet小波變換程序為例說明利用卷積方法實現(xiàn)連續(xù)小波變換的過程(參見程序注釋)。其中，所用morlet小波的定義為

程序如下：

function wcoefs = mymorletcwt(Sig,Scales,fc,fb)

%==================%

% Continuous Wavelet Transform using Morlet function

%======輸入======%

% Sig: 輸入信號

% Scales: 輸入的尺度序列

% fc: morlet小波中心頻率 (默認(rèn)為2)

% fb: morlet小波帶寬參數(shù) (默認(rèn)為2)

%======輸出======%

% wcoefs: morlet小波變換計算結(jié)果

%==================%

if (nargin <= 1),

error('At least 2 parameters required');

end;

if (nargin < 4),

fb = 2;

elseif (nargin < 3),

fc = 2;

end;

wavsupport=8;

% 默認(rèn)morlet小波的支撐區(qū)為[-8,8]

nLevel=length(Scales);

% 尺度的數(shù)目

SigLen = length(Sig);

% 信號的長度

wcoefs = zeros(nLevel, SigLen);?

% 分配計算結(jié)果的存儲單元

for m = 1:nLevel

% 計算各尺度上的小波系數(shù)

a = Scales(m);

% 提取尺度參數(shù)

t = -round(a*wavsupport):1:round(a*wavsupport);

% 在尺度a的伸縮作用下，此時小波函數(shù)的支撐區(qū)會變?yōu)閇-a*wavsup,a*wavsup]，采樣頻率為1Hz

Morl = fliplr((pi*fb)^(-1/2)*exp(i*2*pi*fc*t/a).*exp(-t.^2/(fb*a^2)));

% 計算當(dāng)前尺度下的小波函數(shù)，按小波變換的定義這里需要倒置

temp = conv(Sig,Morl) / sqrt(a);

% 計算信號與當(dāng)前尺度下小波函數(shù)的卷積

d=(length(temp)-SigLen)/2;

% 由于卷積計算所得結(jié)果的長度可能遠(yuǎn)遠(yuǎn)大于原信號，只需提取按原信號的長度提取中間部分的系數(shù)

first = 1+floor(d);

% 區(qū)間的起點

wcoefs(m,:)=temp(first:first+SigLen-1);

end

從以上程序可以看出，基于卷積原理的連續(xù)小波變換實現(xiàn)的關(guān)鍵是求出某一尺度下的小波函數(shù)離散化后的序列。該序列可以通過對該尺度下的小波函數(shù)進(jìn)行采樣求得，采樣區(qū)間為此時小波函數(shù)的支撐區(qū)，采樣頻率可取為1Hz。

注：(1)按我的理解，采樣頻率取得越大，計算結(jié)果也越精確，但我在測試中發(fā)現(xiàn)，采樣頻率取得太高反而會影響分析結(jié)果的精度，在本例中采樣頻率最好取為1Hz。

(2)在小波工具箱的cwt函數(shù)中，某一尺度下的小波函數(shù)采樣序列是直接對母小波采樣序列伸縮而得。

下面利用zhlong給出的例子對mymorletcwt進(jìn)行測試，并與小波工具箱中cwt所得結(jié)果進(jìn)行對比。

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小波進(jìn)行實驗，帶寬參數(shù)為4

fc = 2;

% 中心頻率2Hz

totalscal = 512;

% 所取尺度的數(shù)目

FreqBins = linspace(fmin,fmax,totalscal);

% 將頻率軸在分析范圍內(nèi)等間隔劃分

Scales = fc./ FreqBins;

% 計算相應(yīng)的尺度參數(shù)

wcoefs = mymorletcwt(sig,Scales,fc,fb);

RealFreqBins = FreqBins * SampFreq;

% 尺度所對應(yīng)的實際頻率

%====本帖方法的結(jié)果====%

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結(jié)果==%

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');

測試結(jié)果如下：

本帖方法的結(jié)果

小波工具箱中的cwt結(jié)果

對比兩圖可見，兩種方法的精度相仿。我通過查看cwt的源代碼發(fā)現(xiàn)，兩者的區(qū)別在于小波工具箱中的cwt是按先對小波函數(shù)積分，卷積后再微分的方法來求小波系數(shù)的。至于這樣做的根據(jù)就不得而知了。對于其它小波，也應(yīng)該可以按以上原理實現(xiàn)，特別是那些可以用解析形式表示的小波，如Gaussian小波、Shannon小波、B樣條小波等。

小結(jié)：按卷積原理實現(xiàn)的連續(xù)小波變換的方法簡單直觀，并且個人認(rèn)為其本質(zhì)上是一種矩形數(shù)值積分法。但這種方法的計算精度和速度可能不如其它方法，如更高精度的數(shù)值積分法、調(diào)頻Z變換法、梅林變換法等。在此，偶也希望其它版友能積極發(fā)帖對這幾種方法進(jìn)行討論。

如何利用Matlab進(jìn)行小波分析

?小波分析(wavelet?analysis),?或小波轉(zhuǎn)換(wavelet?transform)是指用有限長或快速衰減的、稱為母小波(mother?wavelet)的振蕩波形來表示信號。該波形被縮放和平移以匹配輸入的信號。

?小波變換的概念是由法國從事石油信號處理的工程師J.Morlet在1974年首先提出的，通過物理的直觀和信號處理的實際經(jīng)驗的需要建立了反演公式，當(dāng)時未能得到數(shù)學(xué)家的認(rèn)可。正如1807年法國的熱學(xué)工程師J.B.J.Fourier提出任一函數(shù)都能展開成三角函數(shù)的無窮級數(shù)的創(chuàng)新概念未能得到著名數(shù)學(xué)家J.L.Lagrange，P.S.Laplace以及A.M.Legendre的認(rèn)可一樣。幸運的是，早在七十年代，A.Calderon表示定理的發(fā)現(xiàn)、Hardy空間的原子分解和無條件基的深入研究為小波變換的誕生做了理論上的準(zhǔn)備，而且J.O.Stromberg還構(gòu)造了歷史上非常類似于當(dāng)前的小波基;1986年著名數(shù)學(xué)家Y.Meyer偶然構(gòu)造出一個真正的小波基，并與S.Mallat合作建立了構(gòu)造小波基的統(tǒng)一方法加多尺度分析之后，小波分析才開始蓬勃發(fā)展起來，其中比利時女?dāng)?shù)學(xué)家I.Daubechies撰寫的《小波十講(Ten Lectures on Wavelets)》對小波的普及起了重要的推動作用。它與Fourier變換、窗口Fourier變換(Gabor變換)相比，這是一個時間和頻率的局域變換，因而能有效的從信號中提取信息，通過伸縮和平移等運算功能對函數(shù)或信號進(jìn)行多尺度細(xì)化分析(MultiscaleAnalysis)，解決了Fourier變換不能解決的許多困難問題，從而小波變化被譽為"數(shù)學(xué)顯微鏡"，它是調(diào)和分析發(fā)展史上里程碑式的進(jìn)展。

本文以《Wavelet transforms and their applications toturbulence》和《A practicalguide to wavelet analysis》兩篇文章的研究為基礎(chǔ)，利用MATLAB，結(jié)合西北地區(qū)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

具體分析：西北地區(qū)的年均ET0存在2-3a的顯著性震蕩周期和6a的準(zhǔn)周期震蕩。對于2-3a的顯著周期震蕩，分別在20世紀(jì)60年代初期到60年代末期和20世紀(jì)80年代末期到21世紀(jì)初期較為顯著，并且處于低值期，2005年以后進(jìn)入減小期；對于6a的準(zhǔn)周期震蕩而言，分別在20世紀(jì)60年代后期到20世紀(jì)70年代前期、20世紀(jì)70年代后期到20世紀(jì)80年代初期和20世紀(jì)80年代中期到21世紀(jì)初期，并且目前處于自2006年以來ET0的低值區(qū)。值得注意的是，在超過28年的分析中，等值線幾乎全為紅色，且沒有閉合的中心，這主要是由于ET0的時間序列僅為56年，超過28年的周期不能明顯地表示出來。

參考文獻(xiàn)

[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.

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。