信号小波分析平面图谱

yzhan79

　　对信号做cwt想得到时间频率图类似伪彩色图一样的平面图谱，不知道怎么操作?

Eminem

问题能说的再详细一点吗 谢谢

猫头鹰先生

请参考以下说明与程序：
此程序实现的是将50HZ和100HZ正弦波分离的程序。并对小波谱分析的实现给出了启发性程序。由于最近太忙，无法解释其中的理论，紧把程序和图形贴出来，供大家讨论和批评指正。

程序如下：
%% 本程序实现小波用作谱分析

clc;clear;

%% 1.正弦波定义
f1=50; % 频率1
f2=100; % 频率2
fs=2*(f1+f2); % 采样频率
Ts=1/fs; % 采样间隔
N=120; % 采样点数

n=1:N;
y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦波混合

figure(1)
plot(y);
title('两个正弦信号')

figure(2)
stem(abs(fft(y)));
title('两信号频谱')

%% 2.小波滤波器谱分析
h=wfilters('db30','l'); % 低通
g=wfilters('db30','h'); % 高通

h=[h,zeros(1,N-length(h))]; % 补零（圆周卷积，且增大分辨率变于观察）
g=[g,zeros(1,N-length(g))]; % 补零（圆周卷积，且增大分辨率变于观察）

figure(3);
stem(abs(fft(h)));
title('低通滤波器图')

figure(4);
stem(abs(fft(g)));
title('高通滤波器图')

%% 3.MALLET分解算法(圆周卷积的快速傅里叶变换实现)
sig1=ifft(fft(y).*fft(h)); % 低通(低频分量)
sig2=ifft(fft(y).*fft(g)); % 高通(高频分量)

figure(5); % 信号图
subplot(2,1,1)
plot(real(sig1));
title('分解信号1')

subplot(2,1,2)
plot(real(sig2));
title('分解信号2')

figure(6); % 频谱图
subplot(2,1,1)
stem(abs(fft(sig1)));
title('分解信号1频谱')

subplot(2,1,2)
stem(abs(fft(sig2)));
title('分解信号2频谱')

%% 4.MALLET重构算法
sig1=dyaddown(sig1); % 2抽取
sig2=dyaddown(sig2); % 2抽取

sig1=dyadup(sig1); % 2插值
sig2=dyadup(sig2); % 2插值
sig1=sig1(1,[1:N]); % 去掉最后一个零
sig2=sig2(1,[1:N]); % 去掉最后一个零

hr=h(end:-1:1); % 重构低通
gr=g(end:-1:1); % 重构高通
hr=circshift(hr',1)'; % 位置调整圆周右移一位
gr=circshift(gr',1)'; % 位置调整圆周右移一位

sig1=ifft(fft(hr).*fft(sig1)); % 低频
sig2=ifft(fft(gr).*fft(sig2)); % 高频
sig=sig1+sig2; % 源信号

%% 5.比较
figure(7);
subplot(2,1,1)
plot(real(sig1));
title('重构低频信号');
subplot(2,1,2)
plot(real(sig2));
title('重构高频信号');

figure(8);
subplot(2,1,1)
stem(abs(fft(sig1)));
title('重构低频信号频谱');
subplot(2,1,2)
stem(abs(fft(sig2)));
title('重构高频信号频谱');

figure(9)
plot(real(sig),'r','linewidth',2);
hold on;
plot(y);
legend('重构信号','原始信号')
title('重构信号与原始信号比较')




顺便说一下，这是分解两个信号的。要分解多个信号，可以采用多层，但注意2抽取的危险。还有一个方法，先为低频2个，高频1个。用一层正反变换。再把低频两个分开，再用一层正反变换。每次都用一层分解和重构。这有些像svm分两类问题，其实也可分三类。

首先，最重要的一点确定采样间隔。50hz，100hz分开它？首先明确子代滤波概念。50hz要用低通，100hz要用高通。因此最佳分割点75hz。因为小波二分频带，所以最高频率150hz。二分后[0,75]，[75，150]。但采样频率fs=2*f=300hz。

其次，什么样的小波。首先最主要的支撑长度，为什么？那就如同说，FIR滤波器为什么高阶效果好一样。所以我选了DB30，支撑长度60。再次，最好线性相位。最后，怎么看滤波器H(N),G(N)频谱。把原滤波器补零。补多少零，分析的信号长度减去滤波其长度。

最后，进行完MALLET分解后，也就是经过低通和高通滤波后。我们的要求达到了。可以通过小波系数的福利叶变换，看频谱。

这和福利叶的滤波器一样的。但小波变换始终是二分的，而福利叶变换是点频的。

关于多层特别小心，因为2抽取会使频谱混叠。所以最好采用些不抽取的算法。
还有小波滤波器是有过渡带的，在75HZ附近的信号频谱是危险的，它会混到两频带中。所以，任何小波分析前，要对信号频谱份量进行分析，任何盲目的分解都是不可行的。

sd

1. cwt函数功能：实现一维连续小波变换的函数。
   
2. cwt函数语法格式： COEFS=cwt(S, SCALES, 'wname')
   
3. COEFS=cwt(S, SCALES, 'wname', 'plot')
   
4. COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE')
   
5. COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM)
   
6. 使用说明：cwt为一维小波变换的函数。
   
7. 
   
8. 格式 COEFS=cwt(S, SCALES, 'wname') 采用'wname'小波，在正、实尺度SCALES下计算向量一维小波系数。
   
9. 格式 COEFS=cwt(S, SCALES, 'wname', 'plot') 除了计算小波系数外，还加以图形显示。
   
10. 格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 计算并画出连续小波变换的系数，并使用PLOTMODE对图形着色。
    
11. 格式 COEFS=cwt(S, SCALES, 'wname', 'plot') 相当于 格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 中的语法 COEFS=cwt(S, SCALES, 'wname', 'absglb')
    
12. 格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM) 能够计算并画出连续小波变换的系数。系数使用PLOTMODE和XLIM进行着色。其中：XLIM=[x1,x2]，并且有如下关系：1<=x1<=x2<=length(S)。
    
13. 
    
14. MODE值 含义
    
15. 'lvl' scale-by-scale着色模式
    
16. 'glb' 考虑所有尺度的着色模式
    
17. 'abslvl'或'lvlabs' 使用系数绝对值的scale-by-scale着色模式
    
18. 'absglb'或'glbabs' 使用系数绝对值并考虑所有尺度的着色模式
    
19. 
    
20. 
    
21. 例子：cwt函数的应用举例
    
22. t=linspace(-1:1:512);
    
23. s=1-abs(t);
    
24. c=cwt(s,1:32,'cgau4');
    
25. c=cwt(s,[64 32 16:-2:2]);
    
26. c=cwt(s,[3 18 12.9 7 1.5],'db2');
    
27. c=cwt(s,1:64,'sym4','abslvl',[100:400]);

复制代码

zhangzy

猫头鹰先生 发表于 2016-6-22 12:52
请参考以下说明与程序：
此程序实现的是将50HZ和100HZ正弦波分离的程序。并对小波谱分析的实现给出了启发 ...

学习了  真不错