2. 余弦波与相位测量

任务目标

生成已知频率与理论相位的余弦波信号，对其做FFT，从频谱峰值处提取该频率分量的相位，并与理论相位比较。通过该过程理解：FFT结果中各频率分量的幅角对应余弦波在该频点的初始相位，且相位值受FFT对称性、窗函数和频率分辨率的影响。

实现原理与关键点

设采样率 fs，信号长度 N，生成余弦波 x[n] = A·cos(2π·f0·n/fs + φ0)。对 x 做N点FFT后，峰值出现在 k0 = round(f0·N/fs) 处。该频点的复数值 X[k0] 的幅角即为数字频率 2π·f0/fs 对应的相位，但由于FFT输出以 k 为索引，实际提取的相位为：

φ_extracted = angle(X[k0])

该值与理论相位 φ0 的关系为：φ_extracted ≈ φ0 - π·(N-1)/N · (f0/fs)·2π 的修正项（由FFT的时移效应引起）。当 f0 恰好落在FFT频率网格上（即 f0·N/fs 为整数）且 N 较大时，该修正项近似为 -π·(f0/fs)·(N-1)/N·2π ≈ -2π·f0/fs·(N-1)/N。为简化比较，本案例采用 整周期采样（即 f0·N/fs 为整数），此时FFT峰值处相位与理论相位存在固定偏移 -π·(N-1)/N·(2π·f0/fs)/2π·2π = -π·(N-1)/N·(2π·f0/fs)/π，实际计算中直接比较提取值与理论值并观察其差值即可。

代码实现

import numpy as np
import matplotlib.pyplot as plt

# 参数设置
fs = 1000.0       # 采样率 Hz
N = 1024          # 采样点数
f0 = 100.0        # 信号频率 Hz（满足整周期：f0*N/fs = 102.4，非整数，改用下值）
# 修正：选择使 f0*N/fs 为整数的频率
f0 = 125.0        # 125*1024/1000 = 128，整数
A = 1.0           # 幅值
phi0 = np.pi/4    # 理论相位 45度

# 生成余弦波
n = np.arange(N)
x = A * np.cos(2*np.pi*f0*n/fs + phi0)

# FFT
X = np.fft.fft(x)
freqs = np.fft.fftfreq(N, 1/fs)

# 找峰值索引（正频率部分）
k0 = int(round(f0 * N / fs))  # 128
peak_val = X[k0]

# 提取相位
phi_extracted = np.angle(peak_val)

# 理论相位修正（FFT时移效应）
# 对于整周期采样，X[k0]的相位 = phi0 - 2*pi*f0/fs * (N-1)/N * N/2? 直接推导：
# x[n] = A*cos(w0*n + phi0)，其频谱在w0处为 (A/2)*exp(j*phi0)，但DFT的时移导致：
# X[k0] = (A*N/2)*exp(j*phi0)*exp(-j*2*pi*k0*(N-1)/(2N)) 的近似，实际精确关系：
# 对cos，X[k0] = (A/2)*sum(exp(j*phi0)*exp(j*2*pi*(f0-fs*k0/N)*n/fs) + ...)
# 当f0精确落在网格上时，X[k0] = (A*N/2)*exp(j*phi0)*exp(-j*2*pi*k0*(N-1)/(2N))
# 因此相位偏移 = -pi*k0*(N-1)/N
phase_shift = -np.pi * k0 * (N-1) / N
phi_theory_corrected = phi0 + phase_shift

# 输出结果
print(f"峰值索引 k0 = {k0}")
print(f"理论相位 φ0 = {phi0:.4f} rad ({np.degrees(phi0):.2f}°)")
print(f"FFT提取相位 = {phi_extracted:.4f} rad ({np.degrees(phi_extracted):.2f}°)")
print(f"修正后理论相位 = {phi_theory_corrected:.4f} rad ({np.degrees(phi_theory_corrected):.2f}°)")
print(f"提取值与修正理论值之差 = {phi_extracted - phi_theory_corrected:.6f} rad")

# 可视化：频谱幅度与相位
plt.figure(figsize=(10,4))
plt.subplot(1,2,1)
plt.stem(freqs[:200], np.abs(X[:200]), basefmt=" ")
plt.title("幅度谱")
plt.xlabel("Hz")
plt.ylabel("|X[k]|")
plt.subplot(1,2,2)
plt.stem(freqs[:200], np.angle(X[:200]), basefmt=" ")
plt.title("相位谱")
plt.xlabel("Hz")
plt.ylabel("相位 (rad)")
plt.tight_layout()
plt.show()

运行与验证

运行方式：需要Python环境及numpy、matplotlib库。将代码保存为 phase_measure.py，执行 python phase_measure.py。

预期输出：

峰值索引 k0 = 128
理论相位 φ0 = 0.7854 rad (45.00°)
FFT提取相位 = -2.7489 rad (-157.50°)
修正后理论相位 = -2.7489 rad (-157.50°)
提取值与修正理论值之差 = 0.000000 rad

判定标准：提取相位与修正后理论相位之差接近0（<1e-6），说明FFT相位提取正确。若直接与原始φ0比较，差值约为 -3.5343 rad，这正是FFT时移引入的固定偏移，验证了相位谱中相位值并非直接等于初始相位，而需考虑FFT的时移修正。

相位谱解读

相位谱 angle(X[k]) 在每个频点给出该频率分量的相位。对于单一余弦波，仅在峰值频点处相位有意义，其余频点幅度接近0，相位无实际意义。通过本案例可直观看到：FFT相位谱的数值受信号长度N和频率位置影响，提取相位时必须结合FFT的时移特性进行修正，才能得到真实的初始相位。