倒频谱的MATLAB实现
看了论坛中很多的内容，学习了不少，发现倒频谱的内容较小。我根据倒频谱的定义和MATLAB，给出了几种倒频谱的计算，程序如下（可在MATLAB直接使用:handshake ）。但是，从获得的图象还是找不出我要找出的频率。例如频率为5，则对应为0.2秒，可是图形在0.2秒出没有任何特征（频率100也是），是不是程序有问题，请各位提出意见。
sf=1000;
nfft=1000;
x=0:1/sf:1;
y=5*cos(2*pi*5*x).^2+3*cos(2*pi*100*x)+randn(size(x));
t=0:1/sf:(nfft-1)/sf;
nn=1:nfft/4;
subplot(3,1,1)
z=rceps(y);
plot(t(nn),abs(z(nn)));
title('z=rceps(y)')
ylabel('幅值');
grid on;
subplot(3,1,2)
yy= real(ifft(log(abs(fft(y)))));
plot(t(nn),abs(yy(nn)));
title('real(ifft(log(abs(fft(x)))))')
xlabel('时间(s)');
ylabel('幅值');
grid on;
subplot(3,1,3)
w=hanning(nfft);
zw=psd(y,nfft,sf,w,nfft/2);
zzw= real(ifft(log(abs(zw))));
plot(t(nn),zzw(nn));
grid on 原帖由 tingxin81 于 2008-1-26 23:54 发表
倒频谱的MATLAB实现
看了论坛中很多的内容，学习了不少，发现倒频谱的内容较小。我根据倒频谱的定义和MATLAB，给出了几种倒频谱的计算，程序如下（可在MATLAB直接使用:handshake ）。但是，从获得的图象还是找不出我要找出的频率。
例如频率为5，则对应为0.2秒，可是图形在0.2秒出没有任何特征（频率100也是），是不是程序有问题，请各位提出意见。
是方法有问题。频谱是反映时间域上的周期性信号，而倒谱是反映频域上的周期性信号。楼主的信号在时间域上是周期性，在频域上能看出；而该信号在频域上不是周期性的，所以在倒谱域上看不到什么明显的信息。楼主的这种加性信号最方便的是能在频域上把它们分辨（见图）。由于5Hz分量平方了，在频域上表现为该频率的2倍上。 感谢，看来是我理解上的问题。
我的例子中的周期函数的频率分别是10Hz（加了个平方项，转换后其实就是10Hz）和100Hz。我通过FFT或PSD（功率谱密度）都可以识别出两个周期函数的频率，我还以为通过倒频谱
rceps也可以分析两个时域周期函数的频率，看来是理解的问题。
麻烦你构造一个实际的例子（主要是时域信号），我想考察我的算法。
[ 本帖最后由 tingxin81 于 2008-1-27 11:26 编辑 ] 本帖最后由 牛小贱 于 2014-7-3 20:43 编辑
为了方便楼主处理，上传了sunday.txt。下载后把后缀改为.wav。读取数据的程序为：
waveFile='sunday.wav';
=wavread(waveFile);
index1=4000;
nfft=1000;
index2=index1+nfft-1;
y=speech(index1:index2);
本帖最后由 牛小贱 于 2014-7-3 20:45 编辑
根据‘songzy41’提供的资源，编辑M文件如下：
clear
clc
close all hidden
format long
waveFile='sunday.wav';
=wavread(waveFile);
index1=4000;
nfft=1000;
index2=index1+nfft-1;
y=speech(index1:index2);
t=0:1/sf:(nfft-1)/sf;
nn=1:nfft/4;
subplot(3,1,1)
z=rceps(y); %MATLAB提供的倒频谱函数
plot(t(nn),abs(z(nn)));
title('z=rceps(y)')
ylabel('幅值');
grid on;
subplot(3,1,2)
%实倒谱定义一：是通过对时域信号的傅里叶变换的幅值求自然对数，然后再做傅里叶逆变换。
yy= real(ifft(log(abs(fft(y)))));
plot(t(nn),abs(yy(nn)));
title('real(ifft(log(abs(fft(x)))))')
xlabel('时间(s)');
ylabel('幅值');
grid on;
subplot(3,1,3)
%实倒谱定义二：是通过对时域信号的功率谱密度的幅值求自然对数，然后再做傅里叶逆变换。
w=hanning(nfft);
zw=psd(y,nfft,sf,w,nfft/2);
zzw= real(ifft(log(abs(zw))));
plot(t(nn),zzw(nn));
grid on
1.得到信号的倒频谱图如下，从第一、二图可以看出，两图谱图一致，因此，断定MATLAB中的倒频谱函数采用的是定义一，即是通过对时域信号的傅里叶变换的幅值求自然对数，然后再做傅里叶逆变换。这样理解对否？
2.是否从你给的信号中：基频成分1/0.005875＝172Hz；1/0.01163＝86Hz。
我发现最近贴附件后都把原有的文字“冲了”，使得必须再输入一次。
楼主的两点理解是对的，因为是处理语音信号，说明该段语音的音调（基频）是172Hz。
而楼主在用psd做倒谱时犯了点小错误：数据长1000点，求得的psd只有501点，只是单边谱。要构成双边谱再求倒谱，也可得到类似图1和2的图。在
zw=psd(y,nfft,sf,w,nfft/2);
后增加一语句：
zw=;
便能得到下图。
[ 本帖最后由 songzy41 于 2008-1-28 11:42 编辑 ] 感谢songzy41，我解决了问题，也认识到了自己的不足。我现在从事发动机振动分析，做了一些信号分析的工作，感觉深度还是不够。
倒谱是反映频域上的周期性信号，我构造时域信号sin(2*pi*10*n*t)其中，n＝1:m。不知道这个算不算频域中的周期信号，他和倍频的关系是不是一回事？呵呵，又提问了哈！ 原帖由 tingxin81 于 2008-1-28 22:46 发表
倒谱是反映频域上的周期性信号，我构造时域信号sin(2*pi*10*n*t)其中，n＝1:m。不知道这个算不算频域中的周期信号，他和倍频的关系是不是一回事？
当把m个信号叠加在一起时，经倒谱分析是能检测出其基频。但倒谱分析主要是用在同态滤波中，可用于某些信号的解褶积。例如提供的语音信号，就能通过倒谱把频谱中的包络和基频频谱分离，如下图所示。
[ 本帖最后由 songzy41 于 2008-1-29 09:35 编辑 ] 好帖子，我学过倒谱后，也一直迷惑，没有应用起来，8楼能给出一个关于倒谱的成功应用实例么？关注中。。。 本帖最后由 牛小贱 于 2014-7-3 20:46 编辑
我构造了一个信号，M文件如下：
clc;
clear;
fs=5000;
nfft=1024;
n=1:nfft/4; %(1,2,3,…，256)
f=0:fs/nfft:fs/2-fs/nfft; %频率向量－横坐标
i=0;
for k=1:6
for t=0:1/fs:1
i=i+1;
x(1,i)=4*sin(2*pi*25*k*t)*sin(2*pi*350*t)+2.5*sin(2*pi*350*t);
end
end
subplot(3,1,1)
y=abs(fft(x,nfft)); %FFT
plot(f(n),y(n))
xlabel('Frequency');
ylabel('|FFT)|')
grid on
subplot(3,1,2)
w=hanning(nfft);
z=psd(x,nfft,fs,w,nfft/2); %PSD
plot(f(n),abs(z(n)));
xlabel('Frequency(Hz)');
ylabel('|psd()|')
grid on
subplot(3,1,3)
t=1000*(0:1/fs:(nfft-1)/fs);%倒频谱的时间向量－横坐标ms
zw=rceps(x); %rceps
plot(t(n),abs(zw(n)));
xlabel('ms')
ylabel('|Cepstrum|');
grid on;得到如下图：图三中1000/40＝25Hz