如何在 FFT 算法里面精确的锁定一个频率.读者问题答疑

import numpy as np
import matplotlib.pyplot as plt
from math import gcd

def demo_case(N, C, title):
    """
    N: 采样点数
    C: 窗口内周期数（相干采样：C 整数）
    f0 = C/N * Fs
    """
    Fs = 1.0  # 只看相对关系即可，设 Fs=1，f0=C/N
    f0 = C / N * Fs
    t = np.arange(N) / Fs
    x = np.sin(2*np.pi*f0*t)

    # 计算相位分布（mod 2π）以及相位点个数
    phase = (2*np.pi*f0*t) % (2*np.pi)
    unique_phases = len(np.unique(np.round(phase, 12)))  # 去重
    g = gcd(N, C)

    # FFT（相干，无泄漏），看峰是否正落 bin
    X = np.fft.rfft(x) / (N/2)  # 简单归一化
    freqs = np.fft.rfftfreq(N, d=1/Fs)
    mag = 20*np.log10(np.abs(X) + 1e-15)

    # 绘图
    fig, axes = plt.subplots(1, 3, figsize=(13, 3.5))

    # (1) 相位分布在单位圆
    theta = np.linspace(0, 2*np.pi, 400)
    axes[0].plot(np.cos(theta), np.sin(theta))  # 单位圆
    axes[0].plot(np.cos(phase), np.sin(phase), '.', ms=6)
    axes[0].set_aspect('equal', 'box')
    axes[0].set_title(f"{title}\nN={N}, C={C}, gcd(N,C)={g}\nunique phases={unique_phases}")
    axes[0].set_xlabel("cos(phase)")
    axes[0].set_ylabel("sin(phase)")
    axes[0].grid(True)

    # (2) 时间域（显示前两周期）
    T0 = 2 / f0  # 两个周期的时长
    idx = t <= T0
    axes[1].plot(t[idx], x[idx], 'o-')
    axes[1].set_title("Time domain (first 2 periods)")
    axes[1].set_xlabel("time")
    axes[1].set_ylabel("amplitude")
    axes[1].grid(True)

    # (3) 频域
    axes[2].plot(freqs, mag, '-')
    axes[2].axvline(f0, ls='--')
    axes[2].set_xlim(0, 4*f0)  # 看 0~4f0 的范围
    axes[2].set_title("FFT magnitude (dB)")
    axes[2].set_xlabel("frequency (bins*Fs)")
    axes[2].set_ylabel("dB")
    axes[2].grid(True)

    plt.tight_layout()
    plt.show()

# 三个对比用例：
# A) gcd(N,C)>1：相位会重复，unique phases 少（例如 N=64, C=16, gcd=16）
demo_case(N=64, C=16, title="Case A: 相位点强重复（gcd>1）")

# B) 互质：gcd(N,C)=1，相位分布最均匀（N=64, C=15）
demo_case(N=64, C=15, title="Case B: 互质（相位均匀展开）")

# C) C 为质数且与 N 互质（N=64, C=17）
demo_case(N=64, C=17, title="Case C: C取质数且互质（更稳）")