关于小波时频图问题

yzhan79

　　一、绘制原理：
　　1.需要用到的小波工具箱中的三个函数cwt()，centfrq()，scal2frq()
　　(1) COEFS = cwt(x,scales,'wname');
　　函数说明：该函数实现连续小波变换，其中x为输入信号，scales为尺度，wname为小波名称。
　　(2)FREQ = centfrq('wname');
　　函数说明：该函数求以wname命名的母小波的中心频率。
　　(3)F = scal2frq(A,'wname',DELTA);
　　函数说明：该函数能将尺度A转换为伪频率F，其中A为尺度，wname为小波名称，DELTA为采样周期。
　　2.尺度与频率之间的关系
　　在尺度与频率之间仅有一个看近似的映射关系。
　　在小波分析中，设a为尺度，fs为采样频率，Fc为小波中心频率，则尺度a对应的伪频率Fa为Fa=Fc*fs/a。显然，根据采样定理，为使小波尺度图的频率范围为(0,fs/2)，尺度范围应为(2*Fc,inf),其中inf表示为无穷大。在实际应用中，只需取尺度足够大即可。
　　3.尺度序列的确定
　　由式Fa=Fc*fs/a可以看出，为使转换后的频率序列是一等差序列，尺度序列必须取为以下形式：
　　c/totalscal, c/(totalscal-1), ...,c/2,c
　　其中，totalscal是对信号进行小波变换时所用尺度序列的长度(通常需要预先设定好)，c为一常数。而尺度c/totalscal所对应的实际频率应为fs/2，于是可得c=2*Fc*totalscal，于是可得到所需的尺度序列。
　　4.时频图的绘制
　　确定了小波基和尺度后，就可以用cwt求小波系数coefs(系数是复数时要取模)，然后用scal2frq将尺度序列转换为实际频率序列f，最后结合时间序列t，用imagesc(t,f,abs(coefs))便能画出小波时频图。
　　二、应用例子
　　下面给出一实际例子来说明小波时频图的绘制。所取仿真信号是由频率分别为50Hz和100Hz的两个正弦分量所合成的信号。

1. 　　% 小波时频分析
   
2. 
   
3. 　　clc;
   
4. 
   
5. 　　clear;
   
6. 
   
7. 　　close all;
   
8. 
   
9. 　　% 原始信号
   
10. 
    
11. 　　fs=1000;
    
12. 
    
13. 　　f1=50;
    
14. 
    
15. 　　f2=100;
    
16. 
    
17. 　　t=0:1/fs:1;
    
18. 
    
19. 　　s=sin(2*pi*f1*t)+sin(2*pi*f2*t); %原始信号
    
20. 
    
21. 　　figure;
    
22. 
    
23. 　　plot(t, s);
    
24. 
    
25. 　　% 连续小波变换
    
26. 
    
27. 　　wavename='cmor3-3';
    
28. 
    
29. 　　totalscal=256;
    
30. 
    
31. 　　Fc=centfrq(wavename); % 小波的中心频率
    
32. 
    
33. 　　c=2*Fc*totalscal;
    
34. 
    
35. 　　scals=c./(1:totalscal);
    
36. 
    
37. 　　f=scal2frq(scals,wavename,1/fs); % 将尺度转换为频率
    
38. 
    
39. 　　coefs=cwt(s,scals,wavename); % 求连续小波系数
    
40. 
    
41. 　　figure;
    
42. 
    
43. 　　imagesc(t,f,abs(coefs));
    
44. 
    
45. 　　set(gca,'YDir','normal');
    
46. 
    
47. 　　colorbar;
    
48. 
    
49. 　　xlabel('时间 t/s');
    
50. 
    
51. 　　ylabel('频率 f/Hz');
    
52. 
    
53. 　　title('小波时频图');
    

复制代码

　　说明：在这个例子中，最好选用复的morlet小波，其它小波的分析效果不好，而且morlet小波的带宽参数和中心频率取得越大，时频图上反映的时频聚集性越好。
　　三、与其他时频分析对比，如短时傅里叶变换时频图

　　高频时小波效果好，因为小波在高频处分辨率可以自动调整，分辨率高。

1. 　　代码：
   
2. 
   
3. 　　% 时频分析
   
4. 
   
5. 　　clc;
   
6. 
   
7. 　　clear;
   
8. 
   
9. 　　close all;
   
10. 
    
11. 　　% 原始信号
    
12. 
    
13. 　　fs=1000;
    
14. 
    
15. 　　f1=50;
    
16. 
    
17. 　　f2=100;
    
18. 
    
19. 　　t=0:1/fs:1;
    
20. 
    
21. 　　s = sin(2*pi*f1*t)+sin(2*pi*f2*t);%+randn(1, length(t));
    
22. 
    
23. 　　% s = chirp(t,20,0.3,300);
    
24. 
    
25. 　　s = chirp(t,20,1,500,'q');
    
26. 
    
27. 　　figure;
    
28. 
    
29. 　　plot(t, s); %原始信号图
    
30. 
    
31. 　　% 连续小波变换时频图
    
32. 
    
33. 　　wavename='cmor3-3';
    
34. 
    
35. 　　totalscal=256;
    
36. 
    
37. 　　Fc=centfrq(wavename); % 小波的中心频率
    
38. 
    
39. 　　c=2*Fc*totalscal;
    
40. 
    
41. 　　scals=c./(1:totalscal);
    
42. 
    
43. 　　f=scal2frq(scals,wavename,1/fs); % 将尺度转换为频率
    
44. 
    
45. 　　coefs=cwt(s,scals,wavename); % 求连续小波系数
    
46. 
    
47. 　　figure;
    
48. 
    
49. 　　imagesc(t,f,abs(coefs));
    
50. 
    
51. 　　set(gca,'YDir','normal')
    
52. 
    
53. 　　colorbar;
    
54. 
    
55. 　　xlabel('时间 t/s');
    
56. 
    
57. 　　ylabel('频率 f/Hz');
    
58. 
    
59. 　　title('小波时频图');
    
60. 
    
61. 　　% 短时傅里叶变换时频图
    
62. 
    
63. 　　figure
    
64. 
    
65. 　　spectrogram(s,256,250,256,fs);
    
66. 
    
67. 　　% 时频分析工具箱里的短时傅里叶变换
    
68. 
    
69. 　　f = 0:fs/2;
    
70. 
    
71. 　　tfr = tfrstft(s');
    
72. 
    
73. 　　tfr = tfr(1:floor(length(s)/2), :);
    
74. 
    
75. 　　figure
    
76. 
    
77. 　　imagesc(t, f, abs(tfr));
    
78. 
    
79. 　　set(gca,'YDir','normal')
    
80. 
    
81. 　　colorbar;
    
82. 
    
83. 　　xlabel('时间 t/s');
    
84. 
    
85. 　　ylabel('频率 f/Hz');
    
86. 
    
87. 　　title('短时傅里叶变换时频图');
    

复制代码

转自：http://blog.sina.com.cn/s/blog_bdfc1ea80102wiis.html

折戟沉沙boy

楼主你好，我有几个疑问希望您能够解答:(1)morlet小波的带宽参数和中心频率取得越大，时频图上反映的时频聚集性越好,这句话理解成复小波的fb和fc越大越好，有没有一个标准或者经常用的值？
(2)totalscal=256；这个尺度长度选择有什么规定，或者我取信号的总长度可以吗？
希望楼主能够解答。谢谢!

zghitboy

折戟沉沙boy 发表于 2016-4-26 10:47
楼主你好，我有几个疑问希望您能够解答:(1)morlet小波的带宽参数和中心频率取得越大，时频图上反映的时频聚 ...

中心频率和带宽

折戟沉沙boy

zghitboy 发表于 2016-4-26 14:08
中心频率和带宽

好的，谢谢！

GuoCL

好资料，谢谢分享。