查找未压缩单通道PCM数据中特定单频率的最优算法

针对单通道PCM中单个指定频率检测的场景，最优选择是Goertzel算法，相比通用FFT有以下优势：

• 仅针对目标频率做计算，运算量仅为同长度FFT的1/N左右，无任何冗余计算，资源消耗极低
• 实现逻辑简单，仅用到基础的乘加运算，没有FFT所需的复数蝶形运算、位反转排序等复杂逻辑，非常适合BASIC、FORTRAN这类无内置信号处理库的语言实现
• 对实值PCM数据适配性好，不需要额外做复数转换预处理

实现前需要先确认3个基础参数，预计算通用系数：

• 采样率Fs：输入PCM数据的采样频率，单位Hz
• 目标检测频率Ft：需要检测的指定单频率，单位Hz
• 检测块长度N：建议取值为N = round(k * Fs / Ft)，其中k为正整数（一般取1~10即可），N越大检测精度越高，但计算量越大，尽量取整可避免频率泄露
• 预计算通用系数：coeff = 2 * cos(2 * π * Ft / Fs)，单次计算可在整个检测块中复用

处理逻辑全为串行逐样本计算，非常容易用BASIC/FORTRAN实现：

1. 初始化3个临时变量：q0 = 0、q1 = 0、q2 = 0
2. 遍历当前检测块的所有N个PCM样本x[n]（n从0到N-1），循环执行以下操作：
   • q0 = coeff * q1 - q2 + x[n]
   • q2 = q1
   • q1 = q0
3. 遍历完成后计算目标频率的能量值：power = q1² + q2² - coeff * q1 * q2
4. 提前根据实际信号幅度设定判定阈值threshold，若power > threshold则判定当前检测块中存在目标频率，否则不存在

如果需要处理连续PCM流，按固定长度N滑动分块重复上述流程即可，块重叠率可按需调整。

FORTRAN 示例片段 #

! 假设输入PCM为16位有符号整数，采样率8000Hz，检测1000Hz频率，块长度80
REAL :: Fs = 8000.0, Ft = 1000.0, PI = 3.1415926535
INTEGER :: N = 80
REAL :: coeff, q0, q1, q2, power, threshold = 1.0e6 ! 阈值根据实际情况调整
INTEGER :: i, x(80) ! x为存储输入PCM样本的数组

! 预计算系数
coeff = 2.0 * COS(2.0 * PI * Ft / Fs)

! 初始化临时变量
q0 = 0.0
q1 = 0.0
q2 = 0.0

! 逐样本迭代
DO i = 1, N
    q0 = coeff * q1 - q2 + REAL(x(i))
    q2 = q1
    q1 = q0
END DO

! 计算能量并判定
power = q1**2 + q2**2 - coeff * q1 * q2
IF (power > threshold) THEN
    PRINT *, "检测到目标频率"
ELSE
    PRINT *, "未检测到目标频率"
END IF

BASIC 示例片段 #

REM 假设输入PCM为16位有符号整数，采样率8000Hz，检测1000Hz频率，块长度80
Fs = 8000!
Ft = 1000!
N = 80
PI = 3.1415926535!
threshold = 1000000! : REM 阈值根据实际情况调整
DIM x(80) AS INTEGER : REM x为存储输入PCM样本的数组

' 预计算系数
coeff = 2! * COS(2! * PI * Ft / Fs)

' 初始化临时变量
q0 = 0!
q1 = 0!
q2 = 0!

' 逐样本迭代
FOR i = 1 TO N
    q0 = coeff * q1 - q2 + x(i)
    q2 = q1
    q1 = q0
NEXT i

' 计算能量并判定
power = q1^2 + q2^2 - coeff * q1 * q2
IF power > threshold THEN
    PRINT "检测到目标频率"
ELSE
    PRINT "未检测到目标频率"
END IF

• 输入的整数PCM样本可以直接代入计算，不需要提前做归一化，阈值根据实际信号的幅度范围调整即可
• 如果需要同时检测少量（≤5个）不同的目标频率，可对每个频率单独跑一遍Goertzel流程，运算量仍然远低于FFT
• 调用内置余弦函数时注意单位是弧度，不要误用角度计算

内容的提问来源于stack exchange，提问作者gmvoeth