回复 #10 songzy41 的帖子
谢谢你的帮助，我对你分析出的信号做个STFT，看一下能看出时频的关系不。[ 本帖最后由 zhlong 于 2007-7-5 21:57 编辑 ] 本帖最后由 wdhd 于 2016-9-7 13:26 编辑
原帖由 zhangnan3509 于 2007-7-5 20:32 发表
呵呵！能不能把小波包处理的效果图和程序传上来，比较一下啊。咱们版里好像小波包的实例不多呀，不吝赐教啊:handshake
希望能见到一个小波包分析的结果
回复 #15 16443 的帖子
把你的短时傅立叶变换的程序和结果图传上来看看 好不好！让大伙都分析一下 我直接用TFTB工具箱中的TFRSTFT函数进行分析的选用了gauss窗，分析的结果就看不出时频的变化 怎样可以把图4.png中频率轴的显示范围缩小（不是比例尺的变化），事其视觉效果放大，也就是使离f＝0很接近的频率能在较远处显示出来。
不知道我说明白了没有 。:@(
回复 #20 16443 的帖子
set(gca,'Ylim',[]) 我用小波分解后，然后去掉近似部分，只留下细节部分。再做STFT分析的结果如下图，大伙看看这个能反映一个怎样的时频关系？ 本帖最后由 wdhd 于 2016-9-7 13:26 编辑原帖由 16443 于 2007-7-5 22:24 发表
我直接用TFTB工具箱中的TFRSTFT函数进行分析的
选用了gauss窗，分析的结果就看不出时频的变化
楼主的图已在帖子http://www.chinavib.com/forum/thread-47661-1-1.html上给出过了，当时我和版主zhlong都给你提出在ftrstft中窗长取得太短。现在你给出采样频率是2500，而主要信号在4Hz附近，该信号一个周期将约有625个样点，而楼主的程序中取窗长为256，则每一次只取4Hz信号的一部分，所以才会出现如楼主的分析图。以下给出2张图。yl3a中给出FFT的结果，上半部是从0~1250的谱图，下半部是只取低频部分，可看到主峰在4Hz附近；yl3b是用TFRSTFT分析的结果，窗长取1024，用海宁窗。程序如下：
load sig;
x=sig;
N=length(x);
X=fft(x);
fs=2500;
nn=1:N;
ff=(0:N/2)/N*fs;
subplot 211; plot(ff,20*log10(abs(X(1:N/2+1)*2/N)));
axis(); grid;
xlabel('Frequency (Hz)'); ylabel('dB');
subplot 212; plot(ff,20*log10(abs(X(1:N/2+1)*2/N)));
axis(); grid;
xlabel('Frequency (Hz)'); ylabel('dB');
nfft=1024;
figure(2);
h=window(nfft-1,'hamming');
=tfrstft(x,1:N,nfft,h);
imagesc(t,f(1:100)*fs,abs(tfr(1:100,:))); axis('xy'); 另外希望楼主说明主要想要解决什么问题，大家给出的意见也已有不少了，但好象都没有抓住楼主需要解决的问题。
还有本帖子是问提取高频信号。信号的频率范围是0~1250，主要信号在4Hz，楼主指的高频信号是从哪儿开始？因为光说高频信号太泛指了，高频和低频都是一个相对概念。
回复 #23 songzy41 的帖子
很好的例子，刚才i运行了一下报错了，h=window(nfft-1,'hamming');是否应该改成=〉h=window(@hamming,nfft-1);回复 #23 songzy41 的帖子
不过songzy41先生 我发现一个问题 好像h=window(nfft-1,'hamming');应该改成h=window('hamming',nfft-1);才能读过去 window函数中用什么样的格式，按我看完全是由 Matlab的版本所决定的。zhangnan3509版主提到的用h=window('hamming',nfft-1)，好象是在Matlab 6.5中，我这里给出Matlab 6.5的help window
WINDOW Generate a standard windowing function
= WINDOW(TYPE,N or BW,MODE,P) calculates a window function
TYPE is one of:
'blackman'
'kaiser' with parameter P (often called beta) (default P=8)
'gaussian' truncated at P std deviations (default P=3)
'hamming'
'hanning'
'harris3' 3-term balckman-harris with 67dB sidelobes
'harris4' 4-term balckman-harris with 92dB sidelobes
'rectangle'
'triangle'
而Matlab 7.1中，也给出help window
WINDOW Window generation.
H=WINDOW(N,NAME,PARAM,PARAM2)
yields a window of length N with a given shape.
N : length of the window
NAME : name of the window shape (default : Hamming)
PARAM: optional parameter
PARAM2 : second optional parameters
我以上的程序为了能读楼主的mat文件在Matlab 7.1中运行，所以用了h=window(nfft-1,'hamming')的格式。而form提到的h=window(@hamming,nfft-1);我也看到过，不知是在哪一种版本中使用的。 本帖最后由 wdhd 于 2016-9-7 13:27 编辑
首先谢谢楼上各位热心的帮助
原帖由 songzy41 于 2007-7-6 09:25 发表
另外希望楼主说明主要想要解决什么问题，大家给出的意见也已有不少了，但好象都没有抓住楼主需要解决的问题。
还有本帖子是问提取高频信号。信号的频率范围是0~1250，主要信号在4Hz，楼主指的高频信号是从哪儿 ...
1。希望楼主说明主要想要解决什么问题：－》我主要是想看看除了基频4Hz以外，在其他频率成分（1000Hz以下），其振幅随着频率的变化趋势，或者在频率随着时间的变化关系，因为现在也不能确定这两种情况到底那个存在。（也就是非线性振动的时频特性和幅频特性）
2。现在振幅跟频率的变化关系可以反映出来了，其结果的分布见22楼，我发的图片。
3。现在还有一个问题是，我的这种看起来比较明显的振幅频率关系图，是经过先小波分解（去掉低频部分，按照songzy41 的意见），再做STFT分析后得到的。我想问问大伙，在用小波分解后可以直接绘制这种时频分布的三维云图嘛 ？
[ 本帖最后由 16443 于 2007-7-7 12:41 编辑 ] 本帖最后由 wdhd 于 2016-9-7 13:27 编辑
原帖由 16443 于 2007-7-7 12:40 发表
1。希望楼主说明主要想要解决什么问题：－》我主要是想看看除了基频4Hz以外，在其他频率成分（1000Hz以下），其振幅随着频率的变化趋势，或者在频率随着时间的变化关系，因为现在也不能确定这两种情况到底那个存在。（也就是非线性振动的时频特性和幅频特性）
2。现在振幅跟频率的变化关系可以反映出来了，其结果的分布见22楼，我发的图片。
3。现在还有一个问题是，我的这种看起来比较明显的振幅频率关系图，是经过先小波分解（去掉低频部分，按照songzy41 的意见），再做STFT分析后得到的。我想问问大伙，在用小波分解后可以直接绘制这种时频分布的三维云图嘛 ？
1，要解决的问题已非常明确了，“除了基频4Hz以外，其他频率成分（1000Hz以下），其振幅随着频率的变化趋势，或者频率随着时间的变化”。
2，22楼的图片似乎不容易看清楚，高的云片往往把低的挡住了，又不能把幅值大的与频率对应关系看清楚。
3，我认为是可以的，把低频部分去掉得到的波形，再作STFT分析。
我在10楼只取了高频单支d6的波形。下程序是按楼主的要求，先小波分解（去掉低频部分）再做STFT分析，观察其它的频率成分（1000Hz以下）。byl11a是去掉低频部分后小波重构得的波形，byl11b对该波形作STFT分析得到的图。可看到去掉4Hz后最强的是25Hz左右，它的幅值也是随时间变化的，而频率也随时间有变化（200Hz和250Hz附近的信号在0.5秒及1秒后才出现）。
load sig;
s=sig;
=wavedec(s,6,'db3');
N=length(s);
T=4e-4;
fs=1/T;
nn=1:N;
ca6=appcoef(c,l,'db3',6);
m=length(ca6);
c=c';
c1=;
s1 =waverec(c1,l,'db3' );
plot((nn-1)*T,s1,'b'); grid;
xlabel('Time'); ylabel('Amplitude');
figure;
nfft=1024;
h=window(nfft-1,'hamming');
=tfrstft(s1',1:N,nfft,h);
imagesc(t,f(1:820)*fs,abs(tfr(1:820,:))); axis('xy');