Function [WT, FreqBins, Scales] = CWT_Morlet(Sig, WinLen, nLevel);
%============================================================%
%  Continuous Wavelet Transform using Morlet function                
%            Sig : 信號                                          
%    WinLen : 小波函數在尺度參數a=1時的長度   （默認為 10）                 
%      nLevel : 頻率軸劃分區間數   （默認為1024）      
% 
%           WT:  返回的小波變換計算結果
%  FreqBins :  返回頻率軸劃分結果（歸一化頻率，最高頻率為0.5） 
%       Scales:   返回與頻率軸劃分值相對應的尺度劃分 （頻率0.5對應的尺度為1）
%============================================================%

if (nargin == 0),
     error('At least 1 parameter required');
end;

if (nargin < 4),
     iShow = 1;
elseif (nargin < 3), 
     nLevel = 1024;
elseif (nargin < 2),
     WinLen = 10;
end;

Sig = hilbert(real(Sig));                     % 計算信號的解析信號
SigLen = length(Sig);                        % 獲取信號的長度
fmax = 0.5;                                         % 設置最高分析頻率 
fmin = 0.005;                                      % 設置最低分析頻率
FreqBins = logspace(log10(fmin),log10(0.5),nLevel);    % 將頻率軸在分析范圍內等對數坐標劃分
Scales = fmax*ones(size(FreqBins))./ FreqBins;             % 計算響應的尺度參數
omg0 = WinLen / 6;                            % 按給定的小波長度計算相應的小波函數中心頻率
WT = zeros(nLevel, SigLen);              % 分配計算結果的存儲單元 

wait = waitbar(0,'Under calculation, please wait...');
for m = 1:nLevel,
    
    waitbar(m/nLevel,wait);
    a = Scales(m);                                   % 提取尺度參數                              
    t = -round(a*WinLen):1:round(a*WinLen);    
    Morl = pi^(-1/4)*exp(i*2*pi*0.5*t/a).*exp(-t.^2/2/(2*omg0*a)^2);   % 計算當前尺度下的小波函數
    temp = conv(Sig,Morl) / sqrt(a);                                                           % 計算信號與小波函數的卷積  
    WT(m,:) = temp(round(a*WinLen)+1:length(temp)-round(a*WinLen));  
    
end; 
close(wait);

WT = WT / WinLen;