脑电信号中的功率谱密度与微分熵

如果你只想要实现方法，可以点击如右链接：这个链接

主要目标

• 熟悉微分熵求解方法
• 思考微分熵与对数能量谱是否真的等价
• 熟悉与反思功率谱密度的求解（仍有未解决之处）

微分熵求解

首先根据bing，了解微分熵的内涵与求解方式。

在脑电信号分析中，DE特征指的是微分熵（Differential Entropy）。微分熵是香农熵在连续变量上的推广形式，它可以用来衡量脑电信号的复杂度

基于参考文献Differential Entropy Feature for EEG-based Vigilance Estimation，得知脑电信号虽然本身不符合正态分布，但其频域信息符合正态分布。因此可对频域信息使用微分熵。如下述公式，i表示第i个频段， $\sigma^2$表示指定频段的方差。

$$
H_i(X) = \frac{1}{2} log(2 \pi e \sigma_i^2)
$$

$$
\hat \sigma_i^2 = \frac{1}{N} \sum_{n=0}^{N-1} (x_i-\hat\mu)^2
$$

求解方式：快速傅里叶变换+方差

$x_i$指的到底是什么？是fft变换后的每个频率下的幅值么。

论文中认为对数能量谱（也可以说是功率谱）等价于微分熵，是假设了某段频率带宽的平均幅度为0，这才能消掉$\hat\mu$。但某段频率的平均幅度一定为零么。论文原话是DC直流分量被去除后，一段信号频段即为0均值了，个人感觉还是有待商榷啊。

For a frequency band of the EEG signals, as a result of zero mean (DC component is filtered)

论文中的下述公式也读不太懂。Pi等于最左边，按我的理解那就是总的功率谱，等于中间的那就是平均功率谱。这三个都等起来那$x_i$和$X_k$到底有啥区别。
$$
\sum_{n=0}^{N-1} |x_i^2|= \frac {1}{N} \sum_{k=0}^{N-1} |X_k^2| = P_i
$$

功率谱密度求解

能量谱，功率谱与功率谱密度

有文章，即下面引用1说功率谱为功率谱密度的简称，但在scipy.signal.periodogram中，参数scaling可以取density或spectrum，分别表示功率谱密度或功率谱，且不同取值的计算结果是有差异的。因此二者似乎并不相同。

scaling : { 'density', 'spectrum' }, optional
Selects between computing the power spectral density ('density') where Pxx has units of V**2/Hz and computing the power spectrum ('spectrum') where Pxx has units of V**2, if x is measured in V and fs is measured in Hz. Defaults to 'density'

首先结合bing与个人理解给出修正后的概念描述：

• 能量谱（energy spectrum）：描述了信号或时间序列的能量如何随频率分布变化1。能量谱是原信号傅立叶变换的平方1能量谱通常用于分析能量有限的信号)4。
• 功率谱（power spectrum）：与能量谱类似，通常针对功率信号而言，描述的是信号的能量在不同频率下的分布，其单位是原始信号单位的平方（在时域中，信号通常表示为某种物理量（如电压、电流、声压等）随时间的变化。当我们将信号从时域转换到频域时，我们实际上是在查看这些物理量在不同频率下的分布。因此，频域信号的单位与时域信号的单位是相同的）。例如，如果原始信号的单位是伏特（V），那么功率谱的单位就是伏特的平方（V²）。
• 功率谱密度(power spectral density)：定义为单位频带内的信号功率1。计算方法为(傅立叶变换的平方)/ (区间长度)

然而，需要注意的是，虽然我们通常将脑电信号视为功率信号，但在进行脑电信号分析时，我们仍然会计算其能量谱。能量谱描述了脑电信号的能量如何随频率分布，这对于理解脑电信号的频域特性是非常有用的。

关于单位的理解

功率谱密度的单位是V**2/Hz。怎么理解呢，首先，傅里叶变换后的平方，此时的单位是V**2。基于定义，单位频带内的信号功率，就用总功率除以这段带宽内的总频率数，进而得到每个频率点的平均功率。/Hz即表示每个频率点。比如对100个频率点进行平均，那么就是用100个原始频率点的总功率除以100个Hz，即$\frac{power\space V^2}{100 Hz}=\frac{power}{100}V^2/Hz$。类似于100个苹果分给10个人，每个人拿到的苹果数为10，记为10个苹果/人。

dB就是$10*log_2(X)$

尚未统一的求解方式

github上求解方式如下。$x$为给定的一段时域上的信号，$w$为窗函数， $N$是序列$x$的总长度。下例中对所有的频谱带宽计算了功率谱密度。
$$
FFTdata = fft(x*w, N)
$$
$$
magFFTdata = abs(FFTdata[0:N//2])
$$
$$
psd = \frac {sum(magFFTdata^2)} {length \space of \space sum}
$$

而在现有库中，求解方式却有所差异。scipy.signal.periodogram函数的源代码实现了基于快速傅里叶变换（FFT）的功率谱密度（PSD）计算。以下是其主要步骤的概述：

1. 预处理：首先，输入的信号会经过一些预处理步骤，包括检查输入的有效性，以及根据需要对信号进行截断或填充。
2. 应用窗函数：然后，会应用一个窗函数到信号上。窗函数的作用是减少信号边缘的突变对FFT结果的影响。scipy.signal.periodogram默认使用的是矩形窗，也就是不做任何改变。
3. 计算FFT：接着，使用scipy.fftpack.fft函数计算信号的快速傅里叶变换。
4. 计算PSD：然后，计算功率谱密度。这是通过取FFT结果的模平方（即实部平方加虚部平方），然后除以采样频率和窗函数能量（窗函数值的平方和，矩形窗的平方和即为原始序列长度或总频率点数）来完成的。
5. 返回结果：最后，函数返回频率值和对应的功率谱密度值。

验证

import numpy as np
from scipy.signal import periodogram

# 生成一个模拟信号
N = 1000
fs = 1000
signal = np.random.randn(1000)

# 计算功率谱
freq, Pxx = periodogram(signal, fs=fs, scaling="density")

fft = np.fft.fft(signal)
psd = np.abs(fft[0:N//2])**2 / fs / (N//2)

import matplotlib.pyplot as plt
plt.plot(Pxx)
plt.plot(psd+0.005)

使用 FFT 获得功率频谱密度估计 - MATLAB & Simulink - MathWorks 中国 中给出了相似的验证。两者相同的是，同样在计算功率谱密度时多除了一个采样频率，为什么？

在scipy.signal.periodogram函数中，除以窗函数能量和采样频率是为了正确地计算功率谱密度。
除以窗函数能量：窗函数是在进行傅里叶变换之前应用到信号上的。窗函数的目的是减少信号边缘的突变对FFT结果的影响。然而，窗函数也会改变信号的能量。为了纠正这一点，我们需要将FFT的结果除以窗函数的能量。
除以采样频率：这是为了将FFT的结果从"每个样本"转换为"每个频率"。FFT的结果是基于所有样本的，但是我们通常关心的是单位频率（例如，每赫兹或每千赫兹）上的功率。因此，我们需要将FFT的结果除以采样频率，以得到功率谱密度。

以上结果为bing上的回答，但并没有很好地说服我。原因有二：

• 其一为窗函数能量。以矩形窗为例，其实根本没有对信号作任何改变，但却除以了$N//2$，此时已经等价于github以及其他文章中提及的功率谱密度的计算方法了。
• 其二为除以采样频率的依据。感觉只是为了转换而转换。除掉窗函数能量后此时已经是每个频率点上的平均功率了，为啥还要再除以采样频率呢？

虽然但是，两种功率谱密度的结果差异仅仅是常数级别的，相差一个采样率。但又想不出来什么哪个应该是对的。

periodogram中density和spectrum的差异

import numpy as np
from scipy.signal import periodogram
import matplotlib.pyplot as plt

# 生成一个模拟信号
N = 10000
fs = 10
signal = np.random.randn(N)

freq, Pxx1 = periodogram(signal, fs=fs, scaling="density")
freq, Pxx2 = periodogram(signal, fs=fs, scaling="spectrum")
plt.plot(Pxx1)
plt.plot(Pxx2*(N//fs)+0.005)

运行上述代码即可发现图形相同。二者差了一个N//fs，原因不明。单位是如何变换成V**2也解释不清楚。