小波分析(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.