声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 40387|回复: 66

[FFT] 倒频谱的MATLAB实现

  [复制链接]
发表于 2008-1-26 23:54 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
看了论坛中很多的内容,学习了不少,发现倒频谱的内容较小。我根据倒频谱的定义和MATLAB,给出了几种倒频谱的计算,程序如下(可在MATLAB直接使用:handshake )。但是,从获得的图象还是找不出我要找出的频率。
例如频率为5,则对应为0.2秒,可是图形在0.2秒出没有任何特征(频率100也是),是不是程序有问题,请各位提出意见。
  1. sf=1000;               
  2. nfft=1000;               
  3. x=0:1/sf:1;
  4. y=5*cos(2*pi*5*x).^2+3*cos(2*pi*100*x)+randn(size(x));  
  5. t=0:1/sf:(nfft-1)/sf;
  6. nn=1:nfft/4;
  7. subplot(3,1,1)
  8. z=rceps(y);
  9. plot(t(nn),abs(z(nn)));
  10. title('z=rceps(y)')
  11. ylabel('幅值');
  12. grid on;
  13. subplot(3,1,2)
  14. yy= real(ifft(log(abs(fft(y)))));
  15. plot(t(nn),abs(yy(nn)));
  16. title('real(ifft(log(abs(fft(x)))))')
  17. xlabel('时间(s)');
  18. ylabel('幅值');
  19. grid on;
  20. subplot(3,1,3)
  21. w=hanning(nfft);      
  22. zw=psd(y,nfft,sf,w,nfft/2);
  23. zzw= real(ifft(log(abs(zw))));
  24. plot(t(nn),zzw(nn));
  25. grid on
复制代码

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2008-1-27 10:36 | 显示全部楼层
原帖由 tingxin81 于 2008-1-26 23:54 发表
倒频谱的MATLAB实现
看了论坛中很多的内容,学习了不少,发现倒频谱的内容较小。我根据倒频谱的定义和MATLAB,给出了几种倒频谱的计算,程序如下(可在MATLAB直接使用:handshake )。但是,从获得的图象还是找不出我要找出的频率。
例如频率为5,则对应为0.2秒,可是图形在0.2秒出没有任何特征(频率100也是),是不是程序有问题,请各位提出意见。

是方法有问题。频谱是反映时间域上的周期性信号,而倒谱是反映频域上的周期性信号。楼主的信号在时间域上是周期性,在频域上能看出;而该信号在频域上不是周期性的,所以在倒谱域上看不到什么明显的信息。楼主的这种加性信号最方便的是能在频域上把它们分辨(见图)。由于5Hz分量平方了,在频域上表现为该频率的2倍上。
tn2b.jpg

点评

赞成: 4.0
赞成: 4
  发表于 2015-3-16 10:46

评分

1

查看全部评分

 楼主| 发表于 2008-1-27 11:20 | 显示全部楼层
感谢,看来是我理解上的问题。
我的例子中的周期函数的频率分别是10Hz(加了个平方项,转换后其实就是10Hz)和100Hz。我通过FFT或PSD(功率谱密度)都可以识别出两个周期函数的频率,我还以为通过倒频谱
rceps也可以分析两个时域周期函数的频率,看来是理解的问题。
麻烦你构造一个实际的例子(主要是时域信号),我想考察我的算法。

[ 本帖最后由 tingxin81 于 2008-1-27 11:26 编辑 ]
发表于 2008-1-27 14:45 | 显示全部楼层
本帖最后由 牛小贱 于 2014-7-3 20:43 编辑

为了方便楼主处理,上传了sunday.txt。下载后把后缀改为.wav。读取数据的程序为:
  1. waveFile='sunday.wav';
  2. [speech, sf, nbits]=wavread(waveFile);
  3. index1=4000;
  4. nfft=1000;
  5. index2=index1+nfft-1;
  6. y=speech(index1:index2);
复制代码

sunday.txt

30.74 KB, 下载次数: 623

点评

赞成: 4.0
赞成: 4
  发表于 2015-3-6 17:08

评分

1

查看全部评分

 楼主| 发表于 2008-1-27 23:21 | 显示全部楼层
本帖最后由 牛小贱 于 2014-7-3 20:45 编辑

根据‘songzy41’提供的资源,编辑M文件如下:
  1. clear
  2. clc
  3. close all hidden
  4. format long
  5. waveFile='sunday.wav';
  6. [speech, sf, nbits]=wavread(waveFile);
  7. index1=4000;
  8. nfft=1000;
  9. index2=index1+nfft-1;
  10. y=speech(index1:index2);
  11. t=0:1/sf:(nfft-1)/sf;
  12. nn=1:nfft/4;
  13. subplot(3,1,1)
  14. z=rceps(y);                      %MATLAB提供的倒频谱函数
  15. plot(t(nn),abs(z(nn)));
  16. title('z=rceps(y)')
  17. ylabel('幅值');
  18. grid on;
  19. subplot(3,1,2)
  20. %实倒谱定义一:是通过对时域信号的傅里叶变换的幅值求自然对数,然后再做傅里叶逆变换。
  21. yy= real(ifft(log(abs(fft(y)))));
  22. plot(t(nn),abs(yy(nn)));
  23. title('real(ifft(log(abs(fft(x)))))')
  24. xlabel('时间(s)');
  25. ylabel('幅值');
  26. grid on;
  27. subplot(3,1,3)
  28. %实倒谱定义二:是通过对时域信号的功率谱密度的幅值求自然对数,然后再做傅里叶逆变换。
  29. w=hanning(nfft);      
  30. zw=psd(y,nfft,sf,w,nfft/2);
  31. zzw= real(ifft(log(abs(zw))));
  32. plot(t(nn),zzw(nn));
  33. grid on
复制代码


1.得到信号的倒频谱图如下,从第一、二图可以看出,两图谱图一致,因此,断定MATLAB中的倒频谱函数采用的是定义一,即是通过对时域信号的傅里叶变换的幅值求自然对数,然后再做傅里叶逆变换。这样理解对否?
2.是否从你给的信号中:基频成分1/0.005875=172Hz;1/0.01163=86Hz。

倒频谱图

倒频谱图
发表于 2008-1-28 09:40 | 显示全部楼层
我发现最近贴附件后都把原有的文字“冲了”,使得必须再输入一次。
楼主的两点理解是对的,因为是处理语音信号,说明该段语音的音调(基频)是172Hz。
而楼主在用psd做倒谱时犯了点小错误:数据长1000点,求得的psd只有501点,只是单边谱。要构成双边谱再求倒谱,也可得到类似图1和2的图。在
zw=psd(y,nfft,sf,w,nfft/2);
后增加一语句:
zw=[zw' zw(end-1:-1:2)'];
便能得到下图。

[ 本帖最后由 songzy41 于 2008-1-28 11:42 编辑 ]
tn11a.jpg

评分

2

查看全部评分

 楼主| 发表于 2008-1-28 22:46 | 显示全部楼层
感谢songzy41,我解决了问题,也认识到了自己的不足。我现在从事发动机振动分析,做了一些信号分析的工作,感觉深度还是不够。
倒谱是反映频域上的周期性信号,我构造时域信号sin(2*pi*10*n*t)  其中,n=1:m。不知道这个算不算频域中的周期信号,他和倍频的关系是不是一回事?呵呵,又提问了哈!
发表于 2008-1-29 09:34 | 显示全部楼层
原帖由 tingxin81 于 2008-1-28 22:46 发表
倒谱是反映频域上的周期性信号,我构造时域信号sin(2*pi*10*n*t)  其中,n=1:m。不知道这个算不算频域中的周期信号,他和倍频的关系是不是一回事?

当把m个信号叠加在一起时,经倒谱分析是能检测出其基频。但倒谱分析主要是用在同态滤波中,可用于某些信号的解褶积。例如提供的语音信号,就能通过倒谱把频谱中的包络和基频频谱分离,如下图所示。

[ 本帖最后由 songzy41 于 2008-1-29 09:35 编辑 ]
tn32a.jpg

评分

2

查看全部评分

发表于 2008-2-1 11:47 | 显示全部楼层
好帖子,我学过倒谱后,也一直迷惑,没有应用起来,8楼能给出一个关于倒谱的成功应用实例么?关注中。。。
 楼主| 发表于 2008-2-3 20:55 | 显示全部楼层
本帖最后由 牛小贱 于 2014-7-3 20:46 编辑

我构造了一个信号,M文件如下:
  1. clc;
  2. clear;
  3. fs=5000;
  4. nfft=1024;
  5. n=1:nfft/4;         %(1,2,3,…,256)
  6. f=0:fs/nfft:fs/2-fs/nfft;   %频率向量-横坐标
  7. i=0;
  8. for k=1:6
  9.     for t=0:1/fs:1
  10.         i=i+1;
  11.         x(1,i)=4*sin(2*pi*25*k*t)*sin(2*pi*350*t)+2.5*sin(2*pi*350*t);
  12.     end
  13. end
  14. subplot(3,1,1)
  15. y=abs(fft(x,nfft));          %FFT
  16. plot(f(n),y(n))
  17. xlabel('Frequency');
  18. ylabel('|FFT)|')
  19. grid on
  20. subplot(3,1,2)
  21. w=hanning(nfft);
  22. z=psd(x,nfft,fs,w,nfft/2);   %PSD
  23. plot(f(n),abs(z(n)));
  24. xlabel('Frequency(Hz)');
  25. ylabel('|psd()|')
  26. grid on
  27. subplot(3,1,3)
  28. t=1000*(0:1/fs:(nfft-1)/fs);  %倒频谱的时间向量-横坐标ms
  29. zw=rceps(x);                %rceps
  30. plot(t(n),abs(zw(n)));
  31. xlabel('ms')
  32. ylabel('|Cepstrum|');
  33. grid on;
复制代码
得到如下图:图三中1000/40=25Hz
untitled.jpg

untitled.fig

13.85 KB, 下载次数: 107

发表于 2008-3-30 16:40 | 显示全部楼层

谢谢!

你是一位热心的前辈,想你致敬!
发表于 2008-11-18 18:26 | 显示全部楼层

1000/40=25Hz

tingxin81!你好!从你的图上好像看不出40ms这个点啊!倒是0.2ms最明显啊!请问这图怎么看?
发表于 2009-4-13 17:21 | 显示全部楼层

why

为什么倒频谱在零点的值非常大啊?尤其是仿真的冲击信号。???
发表于 2009-5-16 17:30 | 显示全部楼层
:victory: 我对此很感兴趣
发表于 2009-5-19 17:48 | 显示全部楼层
学习了,请问倒谱分析有什么用处?与功率谱分析和频谱分析有什么特点?谢谢
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 01:35 , Processed in 0.108625 second(s), 25 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表