小波时频(尺度)图的绘制原理与实现
应广大版友的需要，下面将介绍小波时频（尺度）图的绘制原理，并举例加以说明。－、绘制原理
1.需要用到的小波工具箱中的三个函数
COEFS = cwt(S,SCALES,'wname')
说明：该函数能实现连续小波变换，其中S为输入信号，SCALES为尺度，wname为小波名称。
FREQ = centfrq('wname')
说明：该函数能求出以wname命名的母小波的中心频率。
F = scal2frq(A,'wname',DELTA)
说明：该函数能将尺度转换为实际频率，其中A为尺度，wname为小波名称，DELTA为采样周期。
注：这三个函数还有其它格式，具体可参阅matlab的帮助文档。
2.尺度与频率之间的关系
设a为尺度，fs为采样频率，Fc为小波中心频率，则a对应的实际频率Fa为
Fa＝Fc×fs/a (1)
显然，为使小波尺度图的频率范围为(0,fs/2)，尺度范围应为(2*Fc,inf),其中inf表示为无穷大。在实际应用中，只需取尺度足够大即可。
3.尺度序列的确定
由式(1)可以看出，为使转换后的频率序列是一等差序列，尺度序列必须取为以下形式：
c/totalscal,...,c/(totalscal-1),c/4,c/2,c (2)
其中，totalscal是对信号进行小波变换时所用尺度序列的长度(通常需要预先设定好)，c为一常数。
下面讲讲c的求法。
根据式(1)容易看出，尺度c/totalscal所对应的实际频率应为fs/2，于是可得
c=2×Fc/totalscal (3)
将式(3)代入式(2)便得到了所需的尺度序列。
4.时频图的绘制
确定了小波基和尺度后，就可以用cwt求小波系数coefs（系数是复数时要取模），然后用scal2frq将尺度序列转换为实际频率序列f，
最后结合时间序列t，用imagesc(t,f,abs(coefs))便能画出小波时频图。
注意：直接将尺度序列取为等差序列，例如1:1:64，将只能得到正确的尺度－时间－小波系数图，而无法将其转换为频率－时间－小波系数图。这是因为此时的频率间隔不为常数。
此时，可通过查表的方法将尺度转化为频率或直接修改尺度轴标注。同理，利用本帖所介绍的方法只能得到频率－时间－小波系数图，不能得到正确的尺度－时间－小波系数图。
二、应用例子
下面给出一实际例子来说明小波时频图的绘制。所取仿真信号是由频率分别为100Hz和200Hz的两个正弦分量所合成的信号。
clear;
clc;
fs=1024; %采样频率
f1=100;
f2=200;
t=0:1/fs:1;
s=sin(2*pi*f1*t)+sin(2*pi*f2*t);%两个不同频率正弦信号合成的仿真信号
%%%%%%%%%%%%%%%%%小波时频图绘制%%%%%%%%%%%%%%%%%%
wavename='cmor3-3';
totalscal=256; %尺度序列的长度，即scal的长度
wcf=centfrq(wavename); %小波的中心频率
cparam=2*wcf*totalscal; %为得到合适的尺度所求出的参数
a=totalscal:-1:1;
scal=cparam./a; %得到各个尺度，以使转换得到频率序列为等差序列
coefs=cwt(s,scal,wavename); %得到小波系数
f=scal2frq(scal,wavename,1/fs); %将尺度转换为频率
imagesc(t,f,abs(coefs)); %绘制色谱图
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('小波时频图');
程序运行结果如下：
说明：(1)应用时只须改变wavename和totalscal两个参数即可。
(2)在这个例子中，最好选用复的morlet小波，其它小波的分析效果不好，而且morlet小波的带宽参数和中心频率取得越大，时频图上反映的时频聚集性越好。
[ 本帖最后由 破凰 于 2007-12-17 17:29 编辑 ]
回复 #1 破凰 的帖子
防止灌水和保护帖子的需要，现在将帖子关闭如有其他问题和建议，请专门发帖说明。:handshake回复 楼主 破凰 的帖子
请教你一个一直困扰我的问题'cmor3-3'这里面的贷款系数和中心频率如何确定呢？楼主 破凰 的帖子
楼主，我按照你所给的上述方法用db小波对一个时变信号做时频法分析，但结果和原文出入很大！能不能给解释呢！初来乍到，传文件是给弄错了，结果给弄到我的空间了，烦劳你看看！先谢啦！ 我觉得这一段写的不是很简单明了，绕晕了
改为这样理解，大家看是不是好一些
Fa=Fc*fs/a;
根据实际信号的频率分布范围确定要选择的频率范围fmin,fmax;
但是对于fs采样频率，只有小于或等于fs/2的频率成分，才满足采样定理，所以设置的频率范围fmax<=fs/2
如果想得到频率为等距间隔，那么直接：
f=linspace(fmin,fmax,n);
a=Fc*fs./f;
这样就可以先确定fmin,fmax 来确定a的范围，而且又保证f是等距间隔
我只是觉得第3段写的不易理解，其它都很好
楼主很多帖子我都看过，很佩服啊，
3.尺度序列的确定
由式(1)可以看出，为使转换后的频率序列是一等差序列，尺度序列必须取为以下形式：
c/totalscal,...,c/(totalscal-1),c/4,c/2,c (2)
其中，totalscal是对信号进行小波变换时所用尺度序列的长度(通常需要预先设定好)，c为一常数。
下面讲讲c的求法。
根据式(1)容易看出，尺度c/totalscal所对应的实际频率应为fs/2，于是可得
c=2×Fc/totalscal (3)
将式(3)代入式(2)便得到了所需的尺度序列。
回复 楼主 破凰 的帖子
楼主，你的例程怎么运行不了，我是复制过去的。??? Attempt to execute SCRIPT cwt as a function.
Error in ==> example_cwt at 15
coefs=cwt(s,scal,wavename); %得到小波系数 回复 cuixuepeng 的帖子 你没有小波工具箱
本帖最后由 xinglong-liu 于 2010-8-30 13:47 编辑
回复 chihewanle 如果想得到频率为等距间隔，那么直接：
f=linspace(fmin,fmax,n);
a=Fc*fs./f;
这样就可以先确定fmin,fmax 来确定a的范围，而且又保证f是等距间隔小波的频宽是随中心分析频率变化的，也就是，中心分析频率越高，频宽就越宽。因此，当采用等间隔的分析频率设置时，会出现分析频率的严重重叠，分析效率降低。
c=2×Fc/totalscal 应该是 c=2×Fc×totalscal 很受用，支持，学习 学习了,新手,所以对于此图的实际意义不太了解,觉得没有那种频率-幅值图给人的视觉效果好,可能我还是不太了解小波的缘故吧.得努力了.
能否给出不使用cwt函数的小波变换，s变换、广义s变换的MATLAB程序，非常感谢 回复 破凰 的帖子
弱弱的问一下，为什么最后信号S的时频图上一个时间对应两个频率，两个正弦信号合成后的信号S的频率要怎么求…… 新手，有个问题想请教一下，你的程序画处理的时频图，为什么频率的坐标是倒着的？谢谢了。