yzhan79 发表于 2016-6-2 15:04

信号小波分析平面图谱

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

Eminem 发表于 2016-6-6 14:41

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

猫头鹰先生 发表于 2016-6-22 12:52

请参考以下说明与程序:
此程序实现的是将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=; % 补零(圆周卷积,且增大分辨率变于观察)
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,); % 去掉最后一个零
sig2=sig2(1,); % 去掉最后一个零

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。二分后,。但采样频率fs=2*f=300hz。

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

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

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

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

sd 发表于 2016-6-22 13:12

cwt函数功能:实现一维连续小波变换的函数。
cwt函数语法格式: COEFS=cwt(S, SCALES, 'wname')
COEFS=cwt(S, SCALES, 'wname', 'plot')
COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE')
COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM)
使用说明:cwt为一维小波变换的函数。

格式 COEFS=cwt(S, SCALES, 'wname') 采用'wname'小波,在正、实尺度SCALES下计算向量一维小波系数。
格式 COEFS=cwt(S, SCALES, 'wname', 'plot') 除了计算小波系数外,还加以图形显示。
格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 计算并画出连续小波变换的系数,并使用PLOTMODE对图形着色。
格式 COEFS=cwt(S, SCALES, 'wname', 'plot') 相当于 格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE') 中的语法 COEFS=cwt(S, SCALES, 'wname', 'absglb')
格式 COEFS=cwt(S, SCALES, 'wname', 'PLOTMODE', XLIM) 能够计算并画出连续小波变换的系数。系数使用PLOTMODE和XLIM进行着色。其中:XLIM=,并且有如下关系:1<=x1<=x2<=length(S)。

MODE值 含义
'lvl' scale-by-scale着色模式
'glb' 考虑所有尺度的着色模式
'abslvl'或'lvlabs' 使用系数绝对值的scale-by-scale着色模式
'absglb'或'glbabs' 使用系数绝对值并考虑所有尺度的着色模式


例子:cwt函数的应用举例
t=linspace(-1:1:512);
s=1-abs(t);
c=cwt(s,1:32,'cgau4');
c=cwt(s,);
c=cwt(s,,'db2');
c=cwt(s,1:64,'sym4','abslvl',);

zhangzy 发表于 2016-12-14 12:58

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

学习了真不错
页: [1]
查看完整版本: 信号小波分析平面图谱