1. 背景：

早前接手一个小项目，使用振动传感器监测风机的开关。

这里不对单片机的使用和加速度传感器的数据采集多做说明。

数据准备：振动采样频率设置为400Hz。采样点为512个点。

由离散点FFT可知，实际可以监测的频率范围（400/2=200Hz）。

2. 期望结果

由离散的512个数据点，算出振动频率以及功率。

实际频率F 和数据点横坐标x 成正比关系：

F = x （400/512）=0.78125 x

3. 原始数据点

512个数据点是传感器连续采集得到的。他们是垂直于振动面（Z轴）的加速度值。用加速度值来描述某一瞬间的振动力。

根据加速度传感器芯片量程的配置的不同，原始数据点的数值也会不同。

1. 如下图所示：取一半的点数（256）显示。配置传感器量程为正负16G。传感器水平放置时，Z轴的值表示地球引力G 。图中在2100左右。

2. 用手均匀摇动，显示如下：

下图共256个点，传感器1s钟输出400个点，也就是说图中的3.5个波实际上用了 256/400 =0.64 s 。

手摇的频率在 3.5/0.64 = 5.46875 Hz ,也就是一秒钟手摇了5.5下。（摇的确实很快）

3. 放在电脑主机的侧面，测量主机振动，观察到的原始波形如下：

因为是放在侧面，Z轴受地球引力作用变小（还剩40左右，没有矫正，不过不影响频率的求取）。

4. FFT算法

4.1 基础了解

要对离散的点进行求取频率的操作，首先想到的就是FFT算法，也就是快速傅里叶变换。高等数学，积分变换中有详细的理论解释。网上也有很多的前辈们做出了出色的整理。我这里就不多重复了，给出我当时直接引用的文章链接：

《C》C语言实现FFT算法_杨贵安的博客-CSDN博客_c语言fft

文章里面有封装好的代码，且入参说明的也很详细。

4.2 思路整理

将512个原始数据点传递到FFT函数中，得到新的变化后的数组。数组的个数不变，数据由时域变到了频域。

4.3 代码实现

这里的函数原作者，定义形参时，采用了古老的格式。所以不要惊讶。记得包含头文件"math.h" 。

其中pr[512]数组，存512个原始数据，函数结束后，pr中的数据变为变换后的频域数据模长（用来频谱显示）。

/*******************************************************************\double pr[n]    存放n个采样输入的实部，返回离散傅里叶变换的摸double pi[n]  存放n个采样输入的虚部double fr[n] 返回离散傅里叶变换的n个实部double fi[n]  返回离散傅里叶变换的n个虚部int n 采样点数int k   满足n=2^k
\*******************************************************************/void kfft(pr,pi,n,k,fr,fi)int n,k;double pr[],pi[],fr[],fi[];{ int it,m,is,i,j,nv,l0;double p,q,s,vr,vi,poddr,poddi;for (it=0; it<=n-1; it++)  //将pr[0]和pi[0]循环赋值给fr[]和fi[]{ m=it; is=0;for(i=0; i<=k-1; i++){ j=m/2; is=2*is+(m-2*j); m=j;}fr[it]=pr[is]; fi[it]=pi[is];}pr[0]=1.0; pi[0]=0.0;p=6.283185306/(1.0*n);pr[1]=cos(p); //将w=e^-j2pi/n用欧拉公式表示pi[1]=-sin(p);for (i=2; i<=n-1; i++)  //计算pr[]{  p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);pr[i]=p-q; pi[i]=s-p-q;}for (it=0; it<=n-2; it=it+2)  {vr=fr[it];vi=fi[it];fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];}m=n/2; nv=2;for (l0=k-2; l0>=0; l0--) //蝴蝶操作{ m=m/2; nv=2*nv;for (it=0; it<=(m-1)*nv; it=it+nv)for (j=0; j<=(nv/2)-1; j++){ p=pr[m*j]*fr[it+j+nv/2];q=pi[m*j]*fi[it+j+nv/2];s=pr[m*j]+pi[m*j];s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);poddr=p-q; poddi=s-p-q;fr[it+j+nv/2]=fr[it+j]-poddr;fi[it+j+nv/2]=fi[it+j]-poddi;fr[it+j]=fr[it+j]+poddr;fi[it+j]=fi[it+j]+poddi;}}for (i=0; i<=n-1; i++){ pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);  //幅值计算}return;}

4.4 效果

将传感器放在电脑主机的侧面，用于测频率。电脑主机的硬盘转速是7200转每分钟，对应频率120Hz。

可以看到FFT变换得到的频谱可以提取出120Hz的振动。但低频的噪声影响很严重。这对于后期的数据处理，无疑是增加麻烦。

4.5 优缺点

优点：速度快，1s内可以完成计算。

缺点：随机性变化太大，不利于数据处理。

5. 功率谱

5.1 求取方法

实话实说，我也是看了很多对于频谱信号处理的方法，比较了各种方式的稳定性和易实现程度。选择使用功率谱，就是看到功率谱可以去除掉FFT中的随机性。功率谱的图像更加平稳。

网上前辈们给出了各种求取功率谱的方式：

（1）直接法求功率谱（2）相关函数法（3）相关AR模型法（4） BURG法

相关的知识，大家可以自行baiDu。

我选用的是相关函数法。

5.2 相关函数法

先求取自相关函数，再做频域变换。

5.3 代码实现

1. 在实现代码的过程中发现，需要先对原始值进行预处理。

求出512个原始点的平均值，再对每个原始值减去平均值。这样做的目的是，使得原始值数组分布在零线上下，消除直流分量的影响，重点关注值的变化。

如果不进行该操作，求得的功率谱图像数据整体会很大，频率变化影响的程度在图像上就会很微弱。看起来就是一个毫无波澜对数曲线。

2. 求自相关函数

这里是C语言实现matlab中的xcorr函数。

int xcorr(double *corr, double *x, double *y, int iDataN, int iSyncLength)
{double r =0.0;int i=0, j=0;for (i = 0; i < iDataN- iSyncLength+1; i++){r=0;for(j=0; j < iSyncLength; j++){r+=x[i+j]*y[j];}corr[i]=r;}for (i = iDataN- iSyncLength+1; i < iDataN; i++){r=0;for(j=0; j <iDataN- i; j++){      r+=x[i+j]*y[j];}corr[i]=r;}return 0;}

3. 求功率谱

void PowerSpectrum()
{//求自相关函数xcorr(fr,pr,pr,NUM,NUM);memcpy(pr,fr,NUM*sizeof(double));//快速傅里叶变换kfft(pr,pi,NUM,9,fr,fi);//对数变换{uint16_t i;for(i=0;i<NUM/2;i++){pr[i]=10*log(fabs(pr[i+1]));}}}

5.4 效果

1. 将传感器放在静止的桌面上。可以发现功率谱在没有明显的特征。

功率谱如下：

2. 将硬件传感器放在我的笔记本表面，硬盘转速也是7200转/min（频率是120Hz）。