小波能量计算出现问题

kexin

如果我要分析的信号是非平稳信号，用下面的程序求小波能量时频分布，对么？请指教
y1=y1';
WinLen = 10;
t=t';
[WT, FreqBins, Scales] = CWT_Morlet(y1,WinLen,1024);
FreqBins = FreqBins * fs;
Pyy=abs(WT).^2;
imagesec(t,FreqBins,Pyy); %谱图
colormap jet;
采用普通的Pyy=abs(fft(y1,nfft/2)).^2得出的功率谱图和由上面程序算出的不一样，请问是怎么回事？

songzy41

原帖由 kexin 于 2008-8-12 22:41 发表
小波能量计算出现问题
如果我要分析的信号是非平稳信号，用下面的程序求小波能量时频分布，对么？请指教
y1=y1';
WinLen = 10;
t=t';
[WT, FreqBins, Scales] = CWT_Morlet(y1,WinLen,1024);
FreqBins = FreqBins * fs;
Pyy=abs(WT).^2;
imagesec(t,FreqBins,Pyy); %谱图
colormap jet;
采用普通的Pyy=abs(fft(y1,nfft/2)).^2得出的功率谱图和由上面程序算出的不一样，请问是怎么回事？

在CWT_Morlet函数求出的FreqBins，它不是线性刻度。尽管用imagesec作图时得到的Y轴似乎是线性刻度，实际上问题是imagesec无法做非线性刻度。（不知在MATLAB做图中用什么函数可以做出非线性刻度？）
以下我给出了FreqBins样点和频率对应的图。
又，如果按你昨日的帖子求出最大值的时间和频率，用STFT的方法求出相应最大值的时间和频率，它们两者应该很接近的。

kexin

有点越做越糊涂的感觉，前面是对信号去倾和扣除仪器响应处理，如果对处理后的信号直接用Matlab时频工具箱中的tfrsp(y1), tfrstft(y1), tfrscalo(y1) 三个进行直接分析，得出的能量谱密度图和时间-频率-能量谱图，相差太多了。用小波尺度系数模平方求得的能量，最大值所对应的频率明显偏大。请songzy411推荐一种比较好的可以进行时频分析的方法。

songzy41

原帖由 kexin 于 2008-8-13 10:02 发表
有点越做越糊涂的感觉，前面是对信号去倾和扣除仪器响应处理，如果对处理后的信号直接用Matlab时频工具箱中的tfrsp(y1), tfrstft(y1), tfrscalo(y1) 三个进行直接分析，得出的能量谱密度图和时间-频率-能量谱图，相差太多了。用小波尺度系数模平方求得的能量，最大值所对应的频率明显偏大。请songzy411推荐一种比较好的可以进行时频分析的方法。

我对tfrscalo函数不熟悉，无法评论。但你原来的wavelet_morlet，以及tfrstft和tfrsp得到的结果差不多，所以不知你是怎么调用的。我认为关键是参数的选择。
wavelet_morlet函数是按你原来的程序，得到最大值的时间是3.17秒，对应的频率是2.90662Hz。
tfrstft函数用的参数是：
h=hanning(255);
[BT,tt,FreqBins]=tfrstft(y1,1:n,2048,h);
得到最大值的时间是3.17秒，对应的频率是2.97852Hz。
tfrstft函数用的参数是(窗函数同上)：
[BT,tt,FreqBins]=tfrsp(y1,1:n,2048,h);
得到的结果与ftrstft计算的结果相同，最大值的时间是3.17秒，对应的频率是2.97852Hz。
我认为用tfrsp或tfrstft函数都是很好的时频分析方法。

kexin

您说的对，对于时频分析，窗函数选取及一些参数的正确选取是关键。因为短时傅里叶是固定窗函数，对时频分析有些限制。我倒是用过谱图法进行分析，结果还可以。我试了一下你提供的函数，f=(0:nfft/2)*fs/(nfft);但在用imagesc(t,FreqBins,BT)绘图时又出现问题,如图 ，能量怎么这样分布？是我的频率不对？还是哪里出了问题？而且我求的最大频率值为3.125Hz。能提供一下您求最大频率和时间的程序么？谢谢，盼回复

kexin

更正一下，是imagesc(t,f,BT)。不知你频率是怎么选取的？

songzy41

原帖由 kexin 于 2008-8-13 17:46 发表
对于时频分析，窗函数选取及一些参数的正确选取是关键。因为短时傅里叶是固定窗函数，对时频分析有些限制。我倒是用过谱图法进行分析，结果还可以。我试了一下你提供的函数，f=(0:nfft/2)*fs/(nfft);但在用imagesc(t,FreqBins,BT)绘图时又出现问题,如图布？是我的频率不对？还是哪里出了问题？而且我求的最大频率值为3.125Hz。能提供一下您求最大频率和时间的程序么？谢谢，盼回复
更正一下，是imagesc(t,f,BT)。不知你频率是怎么选取的？

不知道你是否注意到，我们使用的是每帧2048点FFT，而在tfrstft后FreqBins却有4096个，同时BT是2048*1000的数组，即FreqBins与BT不能一一对应，FreqBins中的分辨率要比BT中的更小。所以在显示时我不是直接用FreqBins，作了如下的处理：
[BT,tt,FreqBins]=tfrstft(y1,1:n,2048,h);
FreqBins = FreqBins(1:1024) * fs * 2;
Pyy=abs(BT).^2;
imagesc(t,FreqBins,Pyy(1:1024,:));
axis('xy');
得图如下：

[ 本帖最后由 songzy41 于 2008-8-13 20:09 编辑 ]

kexin

[BT,tt,FreqBins]=tfrstft(y1,1:n,2048,h);
得到最大值的时间是3.17秒，对应的频率是2.97852Hz。
实在是搜不到有好的有关求这种函数最大值时对应的变量的帖子，只能再次拜托songzy411了。也不知道怎么回事，做出来的图像的频率好像比你给出的值要小。:@)

kexin

原帖由 songzy41 于 2008-8-13 20:08 发表

不知道你是否注意到，我们使用的是每帧2048点FFT，而在tfrstft后FreqBins却有4096个，同时BT是2048*1000的数组，即FreqBins与BT不能一一对应，FreqBins中的分辨率要比BT中的更小。所以在显示时我不是直接用FreqBin ...

：我查了一下我们使用的FFT点数是1024吧，而且tfrstft后也是1024个。并未出现翻倍。不知道你用的是不是我前面的程序？而且你把FreqBins = FreqBins(1:1024) * fs*2;为什么要乘以2呢？这样以来，求得的频率比原来不成大2倍。弄不清原理，请讲解一下。谢谢

songzy41

原帖由 kexin 于 2008-8-14 09:30 发表
我查了一下我们使用的FFT点数是1024吧，而且tfrstft后也是1024个。并未出现翻倍。不知道你用的是不是我前面的程序？而且你把FreqBins = FreqBins(1:1024) * fs*2;为什么要乘以2呢？这样以来，求得的频率比原来不成大2倍。弄不清原理，请讲解一下。谢谢

我用的tfrstft的格式如下，主要原来想和wavelet_morlet比较，用的是2048点
[BT,tt,FreqBins]=tfrstft(y1,1:n,2048,h);
我在上一帖子中说到了FreqBins有4096个，而BT是2048*1000的数组。(在Workspace中看数组的大小。)如果用归一化频率来看，FreqBins的分辨率是1/4096=0.000244，而BT的分辨率是1/2048=0.000488，FreqBins的分辨率是BT的一半。所以在BT作图时为了能正确表示出BT的分辨率，又用FreqBins，则可以：
1，FreqBins = FreqBins(1:1024) * fs * 2;
2，FreqBins = FreqBins(1:2:2048) * fs ;
这两个结果是一样的。

songzy41

原帖由 kexin 于 2008-8-14 08:50 发表
[BT,tt,FreqBins]=tfrstft(y1,1:n,2048,h);
得到最大值的时间是3.17秒，对应的频率是2.97852Hz。
实在是搜不到有好的有关求这种函数最大值时对应的变量的帖子，只能再次拜托songzy411了。也不知道怎么回事，做出来的图像的频率好像比你给出的值要小。

我把求出y1后的程序列出来：
%y1=y1';     %不要转置了，tfrstft函数中的x是列向量
h=hanning(255);
[BT,tt,FreqBins]=tfrstft(y1,1:n,2048,h);
FreqBins = FreqBins(1:1024) * fs * 2;
Pyy=abs(BT).^2;
imagesc(t,FreqBins,Pyy(1:1024,:));
axis('xy');

PP=Pyy(:);   %变成一维数组
[Pmax,Ploc]=max(PP);
fprintf('一维数组中最大值的位置: %10d\n',Ploc);
Tm=fix(Ploc/2048)+1;
Fm=mod(Ploc,2048);
fprintf('相应时间位置: %10d   频率位置: %10d\n',Tm,Fm);
fprintf('最大值的时间: %8.5f   频率: %8.5f\n',t(Tm),FreqBins(Fm));

kexin

还有一个问题请教：
我想求能量衰减到最大能量的1/3时，所对应的频率和时间。把程序直接改成这样，怎么结果和最大值返回的结果一样呢？
Pyy=abs(BT).^2;
MM=1/3*Pyy;
PP=MM(:);   %变成一维数组
[Pmax,Ploc]=max(PP);
fprintf('一维数组中最大值的位置: %10d\n',Ploc);
Tm=fix(Ploc/2048)+1;
Fm=mod(Ploc,2048);
fprintf('相应时间位置: %10d   频率位置: %10d\n',Tm,Fm);
fprintf('最大值的时间: %8.5f   频率: %8.5f\n',t(Tm),FreqBins(Fm));