对时序信号做FFT求频谱是现在比较通用的做法之一。然而在科研中,我们常常需要报告统计显著性,如何对频谱做统计检验自然就成为了一个重要的问题。一些研究使用置换检验(之前已经在别的文章中提到过多次),即通过随机打乱时序信号再做FFT,求得白噪声假设下的特定频率的幅值分布(注意不是频谱分布,而是单个频率的幅值概率分布),确定实际数据的频率幅值是否达到显著水平。这种方法不需要考虑数据的正态假设,但仔细讨论一下数据的正态性或许有助于我们针对实验目的设计更好的检验方法。
不过我毕竟不是专业学统计或者信号处理的,如有错误还需指正。
特定频率的幅值的分布
(资料图)
先给我的结论:各个频率的幅值的平方(或者说能量、功率)服从自由度为2的卡方分布。
考虑一个时间序列,它属于一个均值为零的白噪声过程,即。由于实际采样中序列长度必定是有限的,不妨令。序列的总能量为
由于标准正态分布的随机变量的和服从卡方分布,这样就有
现在我们换个角度考虑。假设该序列的能量谱,也就是各个频率的能量大小为。受到采样定理的限制,该能量谱只能给出个独立的频率的能量值。如果在MATLAB中做傅里叶变换的话,我们会看到MATLAB给出的频谱是复数数组。其中,复数的模的平方相当于是能量值,而辐角是初相位。
对于每个特定频率,一个正态白噪声过程将给出一个均匀分布的随机初相位和一个服从的随机能量值。这里可以考虑在复平面坐标上解释这一点,也就是从MATLAB给出的频谱的实部和虚部来考虑。比如,实部就是这样一个形式:,其实就是一个余弦周期序列在某个时间点的取值。我们假设了该序列是正态白噪声,自然我们就会联想到也服从正态分布。类似地,应该也服从正态分布。
当我们假设二者都服从,我们会发现:
又由于白噪声序列的能量谱应该是平的,也就是各个频率有相同能量值,于是
它与一致。不过这里我并不能严格推出实虚部的正态性,同样也不确定方差,文章最后我会给出一些说明。
令。现在我们考虑另一个时序过程。假设在某个频率它有一个不同的参数,类似于,我们有
下面留意一下卡方分布的概率密度函数
我们令,自然地可以将这个函数扩展为
那就有
更简单地讲,当,概率密度函数为
单一频率的能量值检验与能量占比检验
基于以上分布,我们可以给出单一频率的能量值参数的一致最优检验。
考虑我们有个时序数据,其某个频率的能量值为。例如我们的原假设和备择假设是:
由于属于指数族分布,其一致最优检验存在且可由下面的充分统计量给出[1]:
在原假设下,我们有
或者说
确定显著性水平,即可通过该分布判断数据是否达到显著。
然而,更多的一种情况是我们并不确定,但我们希望对某个频率的能量占比做检验。自然地,我们有
考虑对应的随机变量,由于和,又由于服从卡方分布的随机变量之商服从分布,令,我们有
注意到和应该是有单调的关系,也就是当增大时,也会增大。特别当序列是白噪声序列时,由于独立频率有个,因此,而。由于二者的单调关系,对的检验就相当于对的检验。当我们只有一条时序数据时,我们只需要确定原假设,就可以根据该分布确定统计显著性。
但是,当我们有多条时序数据时,我们暂时无法给出具体分布,因为暂时不确定服从分布的随机变量之和的分布表达式。一种方法是通过计算机随机模拟出该分布,不过我们这里也给出另一种近似方法:
如果,那么[2]。
于是当,我们有
如果令为个数据的某频率的能量占比,那么我们就有
我们同样只需要确定原假设,就可以根据该分布确定统计显著性。
频谱形态检验
频谱具有多个频率的信息。有时我们希望检验特定的频谱形态,比如某个频率是峰值点。
考虑原假设
而我们的备择假设是特定频谱形态,其中是某个或某些不确定的参数,比如我们允许峰宽有变化,
令,由于,
于是有概率密度函数:
假设我们有个数据,即有。利用似然比检验,令似然比统计量
令
注意到二者是单调的关系,因此我们可以把作为统计量。我们需要在原假设下生成随机数据,模拟得到该统计量在原假设下的分布。
值得一提的是,我们在确定备择假设时,有可能只能大致确定其形态,换句话说就是知道各个之间的比例,但对其绝对大小把握不好(比如不确定是还是)。我们考虑两种备择假设,满足:
在两种备择假设下我们分别有
二者满足这样形式的关系:
二者也是单调的关系。因此,似乎我们只要把握好备择假设的形态特征即可。
由于所面临的具体问题会不同,对于以上方法的统计功效分析这里就不展示了。
对 (3) 和 (4) 的说明
对于一个服从方差为的正态白噪声过程的时间序列,和中我们假定了某个频率的实部或虚部也都服从方差为的正态分布,这大概是有问题的。注意到各个频率的实部之和相当于序列的第一个点,自然有,这与序列的正态白噪声的方差假设矛盾。后续的推导也应该据此修正。不过巧合的是,对目前所给出的统计检验,我使用不同信噪比的正弦周期序列来模拟,其假阳性率没有问题,统计功效表现也很不错,应该是有某个未知系数在推导过程中可以恰好抵消。这一点还需要后续分析。
参考资料:
[1] https://zhuanlan.zhihu.com/p/102762854 及 https://zhuanlan.zhihu.com/p/87520809
[2] https://en.wikipedia.org/wiki/F-distribution