看到一个倒谱分析的方法，网上查看到一些资料，介绍的挺不错的，分享一下：

这个是东南大学机械工程系的网上课件，略显粗糙，可以将就看

倒谱分析

http://me.seu.edu.cn/jmp/jpkc2006/kczd/kc/chap_10/10_6_4.HTM#####

用倒频谱诊断齿轮故障

http://me.seu.edu.cn/jmp/jpkc2006/kczd/kc/chap_10/10_6_4_3.HTM

分离信息通道对信号的影响

http://me.seu.edu.cn/jmp/jpkc2006/kczd/kc/chap_10/10_6_4_4.HTM

这个是一个新浪博客对上面内容的总结，不知道为什么google里搜不着

http://blog.sina.com.cn/s/blog_4b700c4c01017wzb.html

matlab中的函数为

1) 复倒频谱

y=cceps(x)

2) 实倒频谱

y=rceps(x)

倒谱分析的意义

    对于高速大型旋转机械，其旋转状况是复杂的，尤其当设备出现不对中，轴承或齿轮的缺陷、油膜涡动、磨擦、质量不对称等现象时，则振动更为复杂，用一般频谱分析方法已经难于辩识(识别反映缺陷的频率分量)，而用倒频谱，则会增强识别能力。

    齿轮箱的振动是一种复杂的振动，如齿轮箱中轴的转速为f1 ，轴上齿轮的齿数z ，齿轮箱的振动一定包含有f1，啮合频率z *f1 及其它们各阶谐振。

    频率z*f1的等幅振动在旋转一周期间因轴的振动f1而使振幅有变化，即振幅被调制。因此，在齿轮箱振动谱图上除在频率z*f1处有谱线外，还在(z-1)f1和(z+1)f1处有谱线。考虑被f1的高次谐波调制则出现zf1的边频簇。

http://s7/middle/6163bdeb4ba3eb5fa0b96&690

    如果谱值用对数刻度则该曲线更近似周期波，其“周期”为f1。“周期”为f1的成分在齿轮故障诊断中叫故障函数，与齿轮故障有关需要识别出来。

    一般齿轮箱中部有很多齿轮和转轴、因而有很不同的转轴速度和啮合频率。它们之间的相互调制使得率谱图中包含很多大小和周期都不同的周期成分，在功率谱图混在一起很难分离。

    如果对具有边带信号的功率谱本身再作一次谱分析，则能把边带信号分离出来。因为功率谱图中的周期分量在第二次谱分析的谱图中是离散线段，其高度就反映原功率谱中周期分量的大小。这就是倒频谱分析。

http://s15/middle/6163bdeb4ba3eb609aece&690

    另外，倒谱分析在输入特性与系统特性的分离上也有用途。

实例

    装在燃气轮机与发电机之间的大型齿轮箱，齿轮箱的输入轴转速85s^-1，输出轴转速50s^-1，如图分别表示齿轮箱修理前、后的振动速度功率谱。因中A，B，C三个尖峰对应的是齿轮啮合频率、二次谐波、三次谐波。修理前后的功率谱有一定差别，但是边带无法识别。

http://s14/middle/6163bdeb4ba3eb618a32d&690

    下图分别表示修理前、后的振动速度倒频谱。上图表示边带簇倒频率在11.8ms (85Hz与输入轴转频相同) 处十分明显，说明a 图中的边带，主要是由输入轴的回转误差调制形成的。修理后，11.8ms的分量已大大减小，表明输入轴运转状态已大为改善。

http://s3/middle/6163bdeb4ba3eb62c40c2&690

仿真

    做一个自己想的简单仿真吧，这是一个被5Hz倍频调制的100Hz信号，

http://s15/middle/6163bdeb4ba3eb637d9de&690

    其频谱如下

http://s15/middle/6163bdeb4ba3eb646331e&690

    倒频谱分析如下

http://s15/middle/6163bdeb4ba3eb65663ee&690

     可以发现在t=0.2s处有一个小峰值，对应的就是5Hz，即为那个调制频率。

     代码如下

clc
clear all
close all

ts = 0.001;
Fs = 1/ts;
t = 0:ts:1;
L = length(t);

f0 = 100;
fm = 5;
A = 1;

% 原始信号
y = A*(1+cos(2*pi*fm*t)+cos(2*pi*2*fm*t)+cos(2*pi*3*fm*t)+cos(2*pi*4*fm*t)+cos(2*pi*5*fm*t)).*sin(2*pi*f0*t);
figure
plot(t, y)

% 频谱
NFFT = 2^nextpow2(L);
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
figure
plot(f,2*abs(Y(1:NFFT/2+1)))

% 倒频谱
c = rceps(y);
figure
plot(t,c)