特征值、特征向量与快速傅里叶变换相关技术解析

### 特征值、特征向量与快速傅里叶变换相关技术解析 #### 1. 逆迭代法求特征值和特征向量 逆迭代法的基本思想较为简单。考虑线性系统\((A - \tau I) \cdot y = b\)，其中\(b\)是随机向量，\(\tau\)接近矩阵\(A\)的某个特征值\(\lambda\)，那么解\(y\)会接近对应于\(\lambda\)的特征向量。可以通过迭代来不断逼近真实的特征向量，即每次将\(b\)替换为\(y\)，再求解新的\(y\)。 将\(y\)和\(b\)展开为矩阵\(A\)的特征向量\(x_j\)的线性组合： \(y = \sum_{j} \alpha_j x_j\) \(b = \sum_{j} \beta_j x_j\) 代入线性系统可得： \(\sum_{j} \alpha_j (\lambda_j - \tau) x_j = \sum_{j} \beta_j x_j\) 从而得到： \(\alpha_j = \frac{\beta_j}{\lambda_j - \tau}\) \(y = \sum_{j} \frac{\beta_j x_j}{\lambda_j - \tau}\) 若\(\tau\)接近\(\lambda_n\)，且\(\beta_n\)不是过小，那么\(y\)近似于\(x_n\)。每次迭代会使\(\lambda_j - \tau\)在分母中的幂次增加，因此对于分离良好的特征值，收敛速度很快。 在迭代的第\(k\)阶段，求解方程\((A - \tau_k I) \cdot y = b_k\)，其中\(b_k\)和\(\tau_k\)是当前对特征向量和特征值的猜测。将\(b_k\)归一化，使其满足\(b_k \cdot b_k = 1\)。精确的特征向量和特征值满足\(A \cdot x_n = \lambda_n x_n\)，所以\((A - \tau_k I) \cdot x_n = (\lambda_n - \tau_k) x_n\)。 由于\(y\)是对\(x_n\)的改进近似，将其归一化后得到\(b_{k + 1} = \frac{y}{|y|}\)。通过将改进后的猜测\(y\)代入\((A - \tau_k I) \cdot x_n = (\lambda_n - \tau_k) x_n\)，可以得到特征值的改进估计： \(\tau_{k + 1} = \tau_k + \frac{1}{b_k \cdot y}\) 在实际应用中，逆迭代法的实现可能会比较棘手。首先需要考虑何时使用逆迭代法，大部分计算量在于求解线性系统\((A - \tau_k I) \cdot y = b_k\)。一种可能的策略是先将矩阵\(A\)转换为特殊形式，如对称矩阵的三对角形式或非对称矩阵的 Hessenberg 形式，以便更轻松地求解该线性系统。 逆迭代法的具体实施步骤如下： 1. 提供感兴趣的特征值\(\lambda_n\)的初始值\(\tau_0\)。 2. 选择一个随机归一化向量\(b_0\)作为特征向量\(x_n\)的初始猜测，求解\((A - \tau_0 I) \cdot y = b_0\)。 3. 计算新向量\(y\)相对于\(b_0\)的“增长因子”\(|y|\)，理想情况下该因子应较大。 4. 根据以下情况进行处理： - 如果初始增长因子过小，认为随机向量选择不佳，回到步骤 1 重新选择初始向量。 - 如果\(|b_1 - b_0|\)小于某个容差\(\epsilon\)，可以将其作为停止迭代的准则，最多迭代 5 - 10 次。 - 经过几次迭代后，如果\(|b_{k + 1} - b_k|\)没有快速减小，可以尝试根据\(\tau_{k + 1} = \tau_k + \frac{1}{b_k \cdot y}\)更新特征值。如果\(\tau_{k + 1}\)在机器精度下等于\(\tau_k\)，则停止迭代；否则，使用新的特征值开始新一轮迭代。 逆迭代法可能会遇到两种问题： - 多重或紧密间隔的根：对于对称矩阵，逆迭代法可能只能找到一个特征向量。可以通过扰动\(\tau_0\)的最后几位有效数字，然后重复迭代来获得独立的特征向量。 - 非对称矩阵的缺陷情况：除非初始猜测良好，否则增长因子会很小，且迭代无法改善这种情况。此时可以选择随机初始向量，求解一次\((A - \tau I) \cdot y = b\)，一旦某个向量的增长因子可接受就停止。 对于实矩阵具有复共轭特征值对的情况，需要使用复数运算来求解复特征向量。对于中等大小或更大的非对称矩阵，建议使用包含复特征向量计算的 QR 方法。 #### 2. 快速傅里叶变换（FFT）基础 傅里叶变换是一类重要的计算问题，可用于数据处理和分析。一个物理过程可以在时域（用函数\(h(t)\)表示）或频域（用函数\(H(f)\)表示）进行描述，通过傅里叶变换方程可以在这两个域之间进行转换： \(H(f) = \int_{-\infty}^{\infty} h(t) e^{2\pi i f t} dt\) \(h(t) = \int_{-\infty}^{\infty} H(f) e^{-2\pi i f t} df\) 傅里叶变换具有线性性质，即两个函数之和的变换等于它们各自变换之和，常数乘以函数的变换等于该常数乘以函数的变换。 在时域中，函数\(h(t)\)可能具有一些特殊对称性，这些对称性会在频域中体现出相应的关系，如下表所示： |时域性质|频域性质| | ---- | ---- | | \(h(t)\) 为实函数 | \(H(-f) = [H(f)]^*\) | | \(h(t)\) 为虚函数 | \(H(-f) = -[H(f)]^*\) | | \(h(t)\) 为偶函数 | \(H(-f) = H(f)\) | | \(h(t)\) 为奇函数 | \(H(-f) = -H(f)\) | | \(h(t)\) 为实偶函数 | \(H(f)\) 为实偶函数 | | \(h(t)\) 为实奇函数 | \(H(f)\) 为虚奇函数 | | \(h(t)\) 为虚偶函数 | \(H(f)\) 为虚偶函数 | | \(h(t)\) 为虚奇函数 | \(H(f)\) 为实奇函数 | 此外，傅里叶变换还有一些其他基本性质： - 时间缩放：\(h(at) \Leftrightarrow \frac{1}{|a|} H(\frac{f}{a})\) - 频率缩放：\(\frac{1}{|b|} h(\frac{t}{b}) \Leftrightarrow H(bf)\) - 时间平移：\(h(t - t_0) \Leftrightarrow H(f) e^{2\pi i f t_0}\) - 频率平移：\(h(t) e^{-2\pi i f_0 t} \Leftrightarrow H(f - f_0)\) 对于两个函数\(h(t)\)和\(g(t)\)及其傅里叶变换\(H(f)\)和\(G(f)\)，有两个重要的组合： - 卷积：\(g * h = \int_{-\infty}^{\infty} g(\tau) h(t - \tau) d\tau\)，且\(g * h \Leftrightarrow G(f) H(f)\) - 相关：\(Corr(g, h) = \int_{-\infty}^{\infty} g(\tau + t) h(\tau) d\tau\)，且\(Corr(g, h) \Leftrightarrow G(f) H^*(f)\) 信号的总功率在时域和频域的计算结果相同，这就是 Parseval 定理： \(\int_{-\infty}^{\infty} |h(t)|^2 dt = \int_{-\infty}^{\infty} |H(f)|^2 df\) 通常，人们关心频率区间\([f, f + df]\)内的功率，定义单边功率谱密度（PSD）为： \(P_h(f) = |H(f)|^2 + |H(-f)|^2\)，\(0 \leq f < \infty\) 当\(h(t)\)为实函数时，\(P_h(f) = 2 |H(f)|^2\)。 #### 3. 离散采样数据的傅里叶变换 在实际应用中，函数\(h(t)\)通常是在均匀时间间隔\(\Delta\)下进行采样的，采样值序列为\(h_n = h(n\Delta)\)，\(n = \cdots, -3, -2, -1, 0, 1, 2, 3, \cdots\)。采样间隔\(\Delta\)的倒数称为采样率。 对于任何采样间隔\(\Delta\)，存在一个特殊频率\(f_c\)，称为 Nyquist 临界频率，定义为\(f_c = \frac{1}{2\Delta}\)。Nyquist 临界频率之所以重要，有两个原因： - 好消息：采样定理表明，如果连续函数\(h(t)\)的带宽限制在小于\(f_c\)的频率范围内，即\(H(f) = 0\)对于所有\(|f| \geq f_c\)成立，那么\(h(t)\)完全由其采样值\(h_n\)确定，具体公式为： \(h(t) = \Delta \sum_{n = -\infty}^{\infty} h_n \frac{\sin[2\pi f_c (t - n\Delta)]}{\pi (t - n\Delta)}\) - 坏消息：如果连续函数的带宽不限于小于\(f_c\)，那么采样会导致所有位于频率范围\(-f_c < f < f_c\)之外的功率谱密度被错误地移动到该范围内，这种现象称为混叠。要克服混叠，需要知道信号的自然带宽限制，或者通过模拟滤波来强制设定一个已知的限制，然后以足够快的采样率进行采样，确保每个最高频率周期至少有两个采样点。 下面是离散采样数据傅里叶变换过程的 mermaid 流程图： ```mermaid graph TD; A[开始] --> B[确定采样间隔∆和采样值序列hn]; B --> C[计算Nyquist临界频率fc]; C --> D{带宽是否小于fc}; D -- 是 --> E[根据采样定理确定h(t)]; D -- 否 --> F[处理混叠问题]; F --> G[进行离散傅里叶变换]; E --> G; G --> H[得到离散傅里叶变换结果Hn]; H --> I[结束]; ``` 假设我们有\(N\)个连续的采样值\(h_k = h(t_k)\)，\(t_k = k\Delta\)，\(k = 0, 1, 2, \cdots, N - 1\)，为了估计傅里叶变换\(H(f)\)，我们只在离散频率\(f_n = \frac{n}{N\Delta}\)，\(n = -\frac{N}{2}, \cdots, \frac{N}{2}\)处进行估计。 通过将积分近似为离散和，得到离散傅里叶变换： \(H_n = \sum_{k = 0}^{N - 1} h_k e^{2\pi i k n / N}\) 离散傅里叶变换将\(N\)个复数\(h_k\)映射为\(N\)个复数\(H_n\)，它与时间尺度\(\Delta\)无关。离散傅里叶变换与连续傅里叶变换的关系为\(H(f_n) \approx \Delta H_n\)。 离散傅里叶变换的对称性与连续傅里叶变换几乎相同。离散逆傅里叶变换的公式为： \(h_k = \frac{1}{N} \sum_{n = 0}^{N - 1} H_n e^{-2\pi i k n / N}\) 离散形式的 Parseval 定理为： \(\sum_{k = 0}^{N - 1} |h_k|^2 = \frac{1}{N} \sum_{n = 0}^{N - 1} |H_n|^2\) #### 4. 快速傅里叶变换（FFT）算法 传统计算离散傅里叶变换的方法需要\(N^2\)次复数乘法，而快速傅里叶变换（FFT）算法可以在\(O(N \log_2 N)\)次操作内完成计算，大大提高了计算效率。 FFT 算法的核心是 Danielson - Lanczos 引理，它表明长度为\(N\)的离散傅里叶变换可以重写为两个长度为\(\frac{N}{2}\)的离散傅里叶变换之和，其中一个由原始\(N\)个点中的偶数点组成，另一个由奇数点组成： \(F_k = \sum_{j = 0}^{N - 1} e^{2\pi i j k / N} f_j = \sum_{j = 0}^{N/2 - 1} e^{2\pi i k (2j) / N} f_{2j} + \sum_{j = 0}^{N/2 - 1} e^{2\pi i k (2j + 1) / N} f_{2j + 1} = F_e^k + W^k F_o^k\) 其中\(W = e^{2\pi i / N}\)，\(F_e^k\)和\(F_o^k\)分别是由偶数点和奇数点组成的长度为\(\frac{N}{2}\)的离散傅里叶变换。 FFT 算法的结构分为两个部分： - 第一部分：将数据按位反转排序。这一步只需要交换元素对，不需要额外的存储。 - 第二部分：进行\(\log_2 N\)次迭代，依次计算长度为 2、4、8、\(\cdots\)、\(N\)的变换。在每次迭代中，使用两个嵌套的内循环来实现 Danielson - Lanczos 引理。 以下是实现 FFT 的代码： ```c #include <math.h> #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void four1(float data[], unsigned long nn, int isign) { unsigned long n,mmax,m,j,istep,i; double wtemp,wr,wpr,wpi,wi,theta; float tempr,tempi; n=nn << 1; j=1; for (i=1;i<n;i+=2) { if (j > i) { SWAP(data[j],data[i]); SWAP(data[j+1],data[i+1]); } m=n >> 1; while (m >= 2 && j > m) { j -= m; m >>= 1; } j += m; } mmax=2; while (n > mmax) { istep=mmax << 1; theta=isign*(6.28318530717959/mmax); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; for (m=1;m<mmax;m+=2) { for (i=m;i<=n;i+=istep) { j=i+mmax; tempr=wr*data[j]-wi*data[j+1]; tempi=wr*data[j+1]+wi*data[j]; data[j]=data[i]-tempr; data[j+1]=data[i+1]-tempi; data[i] += tempr; data[i+1] += tempi; } wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } mmax=istep; } } ``` #### 5. 其他 FFT 算法 除了上述的 Cooley - Tukey FFT 算法（按时间抽取），还有 Sande - Tukey FFT 算法（按频率抽取），它先对输入数据进行\(\log_2 N\)次迭代，然后将输出值按位反转排序。对于一些应用，如卷积，可以避免所有的位反转操作，但在实际中，位反转操作只占 FFT 操作计数的一小部分，因此节省的计算时间并不多。 还有一类 FFT 算法将初始数据集长度\(N\)不完全细分到长度为 1 的变换，而是细分到其他小的 2 的幂次，如\(N = 4\)（基 - 4 FFT）或\(N = 8\)（基 - 8 FFT）。这些小变换通过高度优化的代码利用特定小\(N\)的特殊对称性，比简单的 FFT 算法快 20% - 30%。 对于长度\(N\)不是 2 的幂次的数据集，也有相应的 FFT 算法，它们通过类似于 Danielson - Lanczos 引理的关系将初始问题细分为更小的问题，但如果\(N\)的最大质因数较大，这种方法的效果会变差。如果\(N\)是质数，则无法细分，计算复杂度为\(O(N^2)\)。Winograd 傅里叶变换算法是一个例外，它通过对小\(N\)的离散傅里叶变换进行高度优化的编码，并采用新的组合子因子的方法，在某些情况下可以比简单的 FFT 算法快 2 倍，但数据索引更复杂，且不能“原地”执行。 最后，数论变换是一类用于快速卷积的变换，它用模某个大质数\(N + 1\)的整数运算代替浮点运算，用模运算等价物代替 1 的\(N\)次根。虽然它不是真正的傅里叶变换，但性质相似，计算速度可能更快，但应用范围相对较窄，主要用于相关性和卷积计算，因为变换本身不易解释为“频率”谱。 #### 6. 实函数的 FFT 及正弦、余弦变换 在实际应用中，经常需要对实值样本进行 FFT。使用完整的复 FFT 算法处理实数据效率较低，有两种更好的方法： - 同时处理两个实函数：将两个实函数分别作为复输入数组的实部和虚部，然后使用 four1 进行变换，最后根据变换的对称性将结果拆分为两个复数组。以下是实现该方法的代码： ```c void twofft(float data1[], float data2[], float fft1[], float fft2[], unsigned long n) { void four1(float data[], unsigned long nn, int isign); unsigned long nn3,nn2,jj,j; float rep,rem,aip,aim; nn3=1+(nn2=2+n+n); for (j=1,jj=2;j<=n;j++,jj+=2) { fft1[jj-1]=data1[j]; fft1[jj]=data2[j]; } four1(fft1,n,1); fft2[1]=fft1[2]; fft1[2]=fft2[2]=0.0; for (j=3;j<=n+1;j+=2) { rep=0.5*(fft1[j]+fft1[nn2-j]); rem=0.5*(fft1[j]-fft1[nn2-j]); aip=0.5*(fft1[j+1]+fft1[nn3-j]); aim=0.5*(fft1[j+1]-fft1[nn3-j]); fft1[j]=rep; fft1[j+1]=aim; fft1[ ```