倒频谱的MATLAB实现

tingxin81

看了论坛中很多的内容，学习了不少，发现倒频谱的内容较小。我根据倒频谱的定义和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

复制代码

songzy41

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

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

tingxin81

感谢，看来是我理解上的问题。
我的例子中的周期函数的频率分别是10Hz（加了个平方项，转换后其实就是10Hz）和100Hz。我通过FFT或PSD（功率谱密度）都可以识别出两个周期函数的频率，我还以为通过倒频谱
rceps也可以分析两个时域周期函数的频率，看来是理解的问题。
麻烦你构造一个实际的例子（主要是时域信号），我想考察我的算法。

[ 本帖最后由 tingxin81 于 2008-1-27 11:26 编辑 ]

songzy41

为了方便楼主处理，上传了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);
   

复制代码

tingxin81

根据‘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。

songzy41

我发现最近贴附件后都把原有的文字“冲了”，使得必须再输入一次。
楼主的两点理解是对的，因为是处理语音信号，说明该段语音的音调（基频）是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 编辑 ]

tingxin81

感谢songzy41，我解决了问题，也认识到了自己的不足。我现在从事发动机振动分析，做了一些信号分析的工作，感觉深度还是不够。
倒谱是反映频域上的周期性信号，我构造时域信号sin(2*pi*10*n*t)  其中，n＝1:m。不知道这个算不算频域中的周期信号，他和倍频的关系是不是一回事？呵呵，又提问了哈！

songzy41

原帖由 tingxin81 于 2008-1-28 22:46 发表
倒谱是反映频域上的周期性信号，我构造时域信号sin(2*pi*10*n*t)  其中，n＝1:m。不知道这个算不算频域中的周期信号，他和倍频的关系是不是一回事？

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

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

hyl2323

好帖子，我学过倒谱后，也一直迷惑，没有应用起来，8楼能给出一个关于倒谱的成功应用实例么？关注中。。。

tingxin81

我构造了一个信号，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

wangbuliao

你是一位热心的前辈，想你致敬！

liyonglong

tingxin81！你好！从你的图上好像看不出40ms这个点啊！倒是0.2ms最明显啊！请问这图怎么看？

redtu160

为什么倒频谱在零点的值非常大啊？尤其是仿真的冲击信号。？？？

xiaokang0405

:victory: 我对此很感兴趣

dsp242009

学习了，请问倒谱分析有什么用处？与功率谱分析和频谱分析有什么特点？谢谢