频域分析元素
1. 简介
频域分析是信号处理中常用的技术。例如，在心率变异（HRV）分析中，通常根据功率谱分析从RR间隔中提取高频（HF）和低频（LF）成分。作为一名工业工程师，我对这个领域不太熟悉，但我需要这个工具来进行研究。尽管有很多开源软件/库可用，我认为了解正在发生的事情仍然是必要的。因此，我花了一些时间进行频域分析，并在这里做了一个小结。
2. 傅立叶变换
傅里叶变换是频域分析的基础。傅里叶变换将信号从时域映射到频域。有四种类型的傅里叶变换：
- Fourier series:
X(k)=1T∫0Tx(t)e−j2πkt/TdtX(k)=\frac{1}{T} \int_0^T x(t) e^{-j2\pi k t/T} dt - 离散时间傅里叶变换:
X(f)=∑n=−∞∞x(n)e−j2πfnX(f)=\sum_{n=-\infty}^{\infty} x(n) e^{-j2\pi fn} - 连续时间傅里叶变换:
X(f)=∫−∞∞x(t)e−j2πftdtX(f)=\int_{-\infty}^{\infty} x(t) e^{-j2\pi ft}dt - 离散傅里叶变换:
X(k)=∑n=0N−1x(n)e−j2πkn/NX(k)=\sum_{n=0}^{N-1} x(n) e^{-j2\pi kn/N}
每种傅立叶变换都有其逆傅立叶变换，根据频率分量重建原始信号。
以傅立叶级数为例。可以证明任何周期函数 x(t) 可以表示为 cos(2πmt/T) 和 sin(2πnt/T) 的线性组合，其中 T 是 x(t) 的周期，即 x(t+T)=x(t) 对于任何 t 成立，m 和 n 是整数。
x(t)=∑m=0∞amcos(2πmt/T)+∑n=1∞bnsin(2πnt/T) x(t)=\sum_{m=0}^{\infty} a_m \cos(2\pi mt/T) + \sum_{n=1}^{\infty} b_n \sin(2\pi nt/T)
这里ama_m和bnb_n是实数。根据欧拉公式
ejx=cosx+jsinx e^{jx}=\cos x + j \sin x
当 j = \sqrt{-1} 时，我们得到
x(t) = \sum_{k=-\infty}^{\infty} X(k) e^{j2\pi k t/T}
现在 x(t)x(t) 被表示为 ej2πkt/T 的线性组合，其中 X(k)X(k) 是系数（复数），k 是整数。回想一下欧拉公式，ej2πkt/T 有周期 T/k 或频率 k/T。因此 x(t)x(t) 已经被分解成不同频率的组件的组合！可以证明系数 X(k)X(k) 可以计算为
X(k)=\frac{1}{T} \int_0^T x(t) e^{-j2\pi k t/T} dt
这个方程是傅立叶级数。因此，傅立叶级数计算每个频率分量的系数。而前一个方程是反傅立叶级数，根据频率分量的系数重构原始时域信号。
其他类型的傅立叶变换类似。每种傅立叶变换的特征列在下表中。
| | 时域 | 频域 | | | | | --------------------------------- | ---------------- | ----------- | ---------- | ----------- | | 傅立叶级数 | 连续 | 周期性 | 离散 | 非周期性 | | 离散时间傅立叶变换 | 离散 | 非周期性 | 连续 | 周期性 | | 连续时间傅立叶变换 | 连续 | 非周期性 | 连续 | 非周期性 | | 离散傅立叶变换 | 离散 | 周期性 | 离散 | 周期性 |
例如，傅立叶级数将连续周期函数 x(t) 作为输入，并产生离散非周期系数 X(k) 作为输出。表格显示了一些有趣的东西。一个域中的离散总是与另一个域中的周期相关联，而一个域中的连续总是与另一个域中的非周期相关联。
由于数字信号在计算机中以离散值存储，离散时间傅立叶变换（DTFT）和离散傅立叶变换（DFT）更常用。DTFT以离散非周期信号作为输入，并将其频率分量作为连续周期序列输出。DFT以离散周期信号作为输入，并将其频率分量作为离散周期信号输出。如果满足某些条件，例如N=2^n，可以使用快速傅立叶变换（FFT）算法高效计算DFT。
Fourier transform can be considered from the point of view 函数基和核方法, which is the concern of my another article.
还有其他变换，如z变换和拉普拉斯变换。它们与傅里叶变换有非常相似的形式。例如，在z变换中
X(z) = \sum_{n=-\infty}^{\infty} x(n) z^{-n}
这里 zz 是一个复数。让
z=ejw z=e^{jw}
我们得到 DTFT。
3. 能谱和功率谱
能谱被定义为特定频率上的能量
Sx(f) = ∣X(f)∣^2
如果我们对傅立叶级数和连续傅立叶级数在(-∞,∞)范围内进行求和或积分，并对DFT和DTFT在":[-1,1"]范围内进行求和，我们可以得到总能量。
然而，有时信号没有傅里叶变换（连续傅里叶变换或DTFT），即∣X(f)∣|X(f)|趋于无穷大，我们将功率谱定义为Sx(f)S_x (f)除以傅里叶变换的积分区间。
Rx(f) = lim T→∞ 12T|F(x)|2 R_x (f)=limT→∞12T|F(x)|^2
应用傅立叶变换将x(t)在[-T,T]上的频率为f的信号F(x)\mathcal{F} (x)。可以证明Rx(f)R_x (f)可以计算为rx(τ)r_x(\tau)的傅立叶变换。
Rx(f) = F(rx(τ)) R_x (f)=\mathcal{F} (r_x(\tau))
哪里
rx(τ)=limT→∞12T∫−TTx(t+τ)x(t)dt r_x(\tau)=\lim_{T \rightarrow \infty} \frac{1}{2T} \int_{-T}^T x(t+\tau) x(t) dt
4. 静止随机过程
直到现在我们一直在讨论确定性信号，即 x(t) 是 t 的一个特定函数。有时需要考虑随机过程，即 x_t 是一个随机变量。这里我使用下标来区分随机过程和确定性信号。
对于随机过程 x1,x2,⋯ ,xN，如果任意两个变量之间的相关性仅取决于时间差
Cor(xn+k,xn)=Cor(x1+k,x1)=r(k) Cor(x_{n+k}, x_n)=Cor(x_{1+k},x_1)=r(k)
对于任何 nn，期望值是确定的
E(xn)=μ E(x_n)=\mu
随机过程是弱平稳或广义平稳（WSS）的。
对于WSS过程，我们不能直接对原始序列应用傅立叶变换，因为它们是随机数。但我们可以对其自相关函数r(k)r(k)应用傅立叶变换。回顾第3节，
rx(k) = lim┬(T→∞) 1/N ∑┬(0)^(N-1) x(t+k)x(t)dt
是均值为零、方差为1的WSS过程的自相关性（对于其他WSS过程，它将成比例）。得到的值
Rx(f) = F(rk) R_x (f)=\mathcal{F} (r_k)
是功率谱。
5. 功率谱估计
功率谱是通过自相关函数（ACF）的傅里叶变换来计算的。然而，在实践中，该序列通常不是无限的，我们很少知道真实的自相关函数。因此，我们需要估计功率谱。功率谱估计有两类方法：非参数方法和参数方法。
在非参数方法中，ACF直接由数据估计。周期图是一种非参数估计，其中
r^(k)=1N∑n=0N−1−kx(n)x(n+k) \hat{r}(k)=\frac{1}{N} \sum _{n=0}^{N-1-k} x(n)x(n+k)
这里 NN 是信号的长度。它等同于无限级数 y(n)y(n) 的 ACF，其中 y(n)=x(n)y(n)=x(n) 当 n=0,1,⋯ ,N−1n=0,1,\cdots,N-1 且 y(n)=0y(n)=0 对于其他 nn。
在参数方法中，首先拟合一个模型，例如自回归移动平均（ARMA）模型。然后可以计算自相关函数（ACF）。
功率谱分析的限制在于过程应该是平稳的。然而，在许多情况下，这个假设并不成立。此外，功率谱分析需要一长串数据才能获得良好的估计。因此，在许多情况下，状态空间模型可能是一个更好的选择。
获取更多信息
关于 ARMA 模型:
- Cryer, J. D., & Chan, K. S. (2008). 时间序列分析：在R中的应用. Springer.
关于傅立叶变换：
关于功率谱：
- Hayes, M. H. (2009). 统计数字信号处理与建模. John Wiley & Sons.