论坛 › 信号处理方法 › 连续小波变换实现方法的小结与程序详解

破凰 发表于 2007-12-23 19:13

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

在帖子“给大家分享我自己编的程序-连续小波变换” 中，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();
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();
ylabel('Frequency / Hz');
xlabel('Time / sec');测试结果如下：
本帖方法的结果

小波工具箱中的cwt结果


   对比两图可见，两种方法的精度相仿。我通过查看cwt的源代码发现，两者的区别在于小波工具箱中的cwt是按先对小波函数积分，卷积后再微分的方法来求小波系数的。至于这样做的根据就不得而知了。对于其它小波，也应该可以按以上原理实现，特别是那些可以用解析形式表示的小波，如Gaussian小波、Shannon小波、B样条小波等。

   小结：按卷积原理实现的连续小波变换的方法简单直观，并且个人认为其本质上是一种矩形数值积分法。但这种方法的计算精度和速度可能不如其它方法，如更高精度的数值积分法、调频Z变换法、梅林变换法等。在此，偶也希望其它版友能积极发帖对这几种方法进行讨论。

[ 本帖最后由 破凰 于 2007-12-23 19:47 编辑 ]

zhangnan3509 发表于 2007-12-23 19:18

回复 #1 破凰 的帖子

保护原创帖，所以将该帖关闭，如果有对其中的问题有兴趣讨论的版友请重新开帖。谢谢

oww 发表于 2011-2-13 21:56

thank you，这个帖子太好了，赞一个

ljj722 发表于 2012-3-16 12:58

构造小波基函数时，采样频率取1不合理吧？用wavemenu命令打开小波工具箱主菜单，进入Complex Continuous Wavelet 1-D，里面有个Sampling Period，不知道采样周期的选取影响什么。另外归一化频率貌似等于信号实际频率除以奈奎斯特频率，也就信号实际频率除以采样率的一半。

ljj722 发表于 2012-3-16 13:01

楼主的程序我试了下，有问题，得到的系数虚部为0了！
下面是我用的测试程序：
clc;clear all;close all;format long;

Fs=100*10^3;                        %采样率
dt=1/Fs;                            %采样时间周期
f0=4*10^3;                        %信号载频

t=0:dt:dt*1000;                     %时间
g=1.333-cos(0.2*2*pi*f0*t+(pi/2));
y1=g.*sin(2*pi*f0*t);               %仿真信号

%方法一：利用mymorletcwt函数文件，尺度scale=1
fb=0.00001^2;
fc=f0;
temp1=mymorletcwt(y1,1,fc,fb);
figure(1);subplot(2,2,1);plot(t,temp1,'b');

%方法二：利用cmorwavf构建Morlet复值小波，尺度scale=1
scale=1;
fb=0.00001^2;
fc=f0;
lb=-8;ub=8;
t2=-8:1:8; %采样率取1
n=length(t2);
=cmorwavf(lb,ub,n,fb,fc); %生成一个Morlet复值小波，n+1个点
figure(2);subplot(211);plot(x,real(psi));title('Complex Morlet wavelet cmor');xlabel('Real part');grid on;
subplot(212);plot(x,imag(psi));xlabel('Imaginary part');grid on;

wcoef=conv(y1,fliplr(psi))/sqrt(scale);         %计算信号与尺度1下小波函数的卷积
d=(length(wcoef)-length(y1))/2;               %卷积计算所得结果的长度大于原信号，根据原信号的长度只提取中间部分的系数
first=1+floor(d);                               %区间的起点
temp2=wcoef(1,first:first+length(y1)-1);
figure(1);subplot(2,2,2);plot(t,temp2);


%方法三：利用matlab小波工具箱中的cwt命令，小波基函数取cmor1-1，fb=1，fc=1
fb=1;
fc=1;
sca=fc/f0/dt;   %尺度
temp3=abs(cwt(y1,sca,'cmor1-1'))*2/sqrt(sca);
subplot(2,2,4);plot(t,g,'g',t,temp3,'b');

mingmingtree 发表于 2012-4-9 15:19

回复 7 # ljj722 的帖子

高手啊！！
请问您的课题也是用小波处理光信号么？
如果可以，希望能向您请教呢！:@)

西门牛 发表于 2013-3-26 21:14

{:{03}:}

和丹 发表于 2013-7-15 11:05

非常有用，学习了，谢谢你{:3_47:}

赵强生 发表于 2013-8-18 17:00

赞楼主，学习了

zkj_yy 发表于 2015-6-24 11:49

点睡的冯绍峰更送发箍是

563524748 发表于 2015-7-7 20:11

正是需要的 感谢楼主的分享！！！

capuchino 发表于 2015-9-2 07:29

非常好，谢谢分享

fucheng 发表于 2015-9-19 11:13

{:{39}:}此贴必火

1713573225 发表于 2016-1-16 21:15

谢谢分享资料

shanzishanxin 发表于 2016-2-7 10:46

学习中，感谢~

查看完整版本: 连续小波变换实现方法的小结与程序详解