原创 连续小波变换实现方法的小结与程序详解

复制内容到剪贴板 代码: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

复制内容到剪贴板 代码: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');