随机振动信号的特征——PSD（功率谱密度）

理论知识

1、背景
信号可分为确定性信号和随机信号。确定性信号是每个时间点上的值可以用某个数学表达式或图标唯一地确定的信号；而随机信号（random signal)，幅度未可预知但又服从一定统计特性的信号，又称不确定信号（百度百科的解释）。随机信号是普遍存在的，也不能用一个确切的数学公式来描述，因为也不能准确进行预测的信号。正是因为随机信号是随机的，所以只能用统计的方法进行描述，在一定的准确性或可信性范围内用统计学规律去表征随机信号的特性。
如果随机信号的概率特性不随时间变化而变化，则成为平稳随机信号。

2、统计特征量

以上这些统计学特征用于分析信号的特征，并计算出PSD，功能如下：
数学期望值，描述随机信号的平均值。
方差值，描述随机信号幅度变化的强度。
概率密度函数，是描述信号振幅数值的概率。
相关函数，描述随机信号的每两个具有一定时间间隔的幅度值之间的联系程度的数值，它是时间间隔的一个函数。
功率谱密度，描述随机信号在平均意义上的功率谱特性。

3、功率谱密度
由中心极限定理：大量相互独立的随机变量，其平均值正态分布。可知平稳随机信号是趋于正态分布的，所以可以用统计学上特征量来描述随机信号。平稳随机信号的数学期望基本为零，数学期望为0时，方差等于均方值，也就是平均功率。为了描述平均功率，需要使用频谱分析。对于一个随机信号而言，时域信息是杂乱无章的，唯一的确定性信息但是在统计意义下得到的，即幅值呈正态分布，均方值也就是平均功率是固定的。

根据帕塞瓦尔定理，信号在时域的总功率等于在频域的总功率，所以只需要对时域信号进行傅里叶变换即可转换到频域，得到频域分布，只需要将频域分布平方积分就可以得到功率谱密度（PSD）。
但是因为时域的信号是随机的，无法用数学表达式描述，随机信号的不满足傅里叶变换绝对值可积的条件，严格意义傅里叶变换不存在，也就无法通过傅里叶变换将时域转成频域。（傅里叶理论认为任何复杂的波形都可以分解成不同正弦波，而不同的正弦波也能叠加成复杂波形。）

为了解决这个问题，引入自相关函数，它反映了随机信号本身在不同时刻的相互关系，再直白一点：把一个信号平移一段距离，跟原来有多相似。自相关函数将信号的蕴含的周期信号识别出来，并将相位信息去掉（相位不影响平均功率）。

维纳-辛钦定理：功率谱密度函数与自相关是傅里叶变换对。即一个信号的功率密度谱，就是其自相关函数的傅里叶变换。

总结：随机信号的幅值是满足正态分布的，用它的自相关函数求均方值即得到功率谱密度（PSD），一个随机的信号（随机振动）用PSD来描述其特征。

理论实践
传感器采集的样本为离散的点，没必要把离散点还原为原始信号，直接将离散的点用普估计方法求PSD。普估计方法由周期图法和welch方法，周期图法属于有偏估计，误差很大，所以采用welch方法，python或matlab有welch方法的库函数，可直接调用求得PSD。

python代码示例

#该加速度传感器的数据是三轴的X/Y/Z，默认Z轴为1g。
#返回四个值分别为频率和三个轴的PSD值，可用于绘图
def get_psd(payload):
	SAMPLES_NUM = 1024
    TicksInSamplingTimePeriod = int(0.5 + (0.0025 / (1 / 32768)))
    Fs_Acc = 32768 / TicksInSamplingTimePeriod  # 400 Hz
    NFFT_Acc = 512  # SAMPLES_NUM_PER_PIN    # 512
    Sensitivity_Acc = 1 / 1024  # LSB = 1/1024 [g] for the +/-2g range
    window = hamming(M=NFFT_Acc)

    samples = np.array(struct.unpack("<{}h".format(SAMPLES_NUM * 3), payload))#从payload中解析出来数据sample
    samples = np.array(samples) * Sensitivity_Acc
    acc_x_units_g, acc_y_units_g, acc_z_units_g = samples[0::3], samples[1::3], samples[2::3]
    (freq_Acc, Pxxf_Acc) = welch(x=acc_x_units_g, fs=Fs_Acc, window=window, noverlap=NFFT_Acc / 2,
                                 nfft=NFFT_Acc, return_onesided=True, detrend=False)
    Pxxf_AccX_dBg = 10. * np.log10(abs(Pxxf_Acc) + epsilon)
    (freq_Acc, Pxxf_Acc) = welch(x=acc_y_units_g, fs=Fs_Acc, window=window, noverlap=NFFT_Acc / 2,
                                 nfft=NFFT_Acc, return_onesided=True, detrend=False)
    Pxxf_AccY_dBg = 10. * np.log10(abs(Pxxf_Acc) + epsilon)
    (freq_Acc, Pxxf_Acc) = welch(x=acc_z_units_g, fs=Fs_Acc, window=window, noverlap=NFFT_Acc / 2,
                                 nfft=NFFT_Acc, return_onesided=True, detrend=False)
    Pxxf_AccZ_dBg = 10. * np.log10(abs(Pxxf_Acc) + epsilon)
    return freq_Acc, Pxxf_AccX_dBg, Pxxf_AccY_dBg, Pxxf_AccZ_dBg