Python与Swift功率谱实现结果差异问题求助

问题背景 #

我有一个实信号HRV(t)，需要计算其功率谱。分别用Python和Swift实现代码后，两者输出结果差异极大。用单频正弦波测试时，Python输出严格的delta函数，Swift却并非如此。

疑问点 #

• 是否存在频谱泄漏？
• 两种实现是否使用不同基函数？
• 数据传入Swift的方式是否有误？

实现代码 #

Python代码 #

import numpy as np

def power_spectrum(signal, fs=30):
    """
    Compute the power spectrum of a 1D real signal.

    Parameters:
    signal (array-like): Input signal
    fs (float): Sampling frequency (default is 30)

    Returns:
    f (numpy array): Frequencies corresponding to the power spectrum
    power (numpy array): Normalized power spectrum
    """
    # Length of the signal
    N = len(signal)

    # Frequency array
    f = np.arange(0, N // 2 + 1) * fs / N

    # Compute the FFT of the signal
    fft_result = np.fft.fft(signal)

    # Compute the power spectrum
    power = (1 / N) * np.abs(fft_result[:N // 2 + 1])**2

    # Normalize the power spectrum
    power = power / np.max(power)

    return f, power

Swift代码 #

public func powerSpectrum(_ input: [Float]) -> [Float] {
    // apply hamming window
    // window has slightly different precision than scipy.window.hamming (float vs double)
    var win = [Float](repeating: 0, count: input.count)
    vDSP_hamm_window(&win, vDSP_Length(input.count), 0)
    var input_windowed = [Float](repeating: 0, count: input.count)
    vDSP_vmul(win, 1, input, 1, &input_windowed, 1, vDSP_Length(win.count))
    symmetricZeroPad(signal: &input_windowed, length: 2048)
    var real = [Float](input_windowed)
    var imaginary = [Float](repeating: 0.0, count: input_windowed.count)
    var splitComplex = DSPSplitComplex(realp: &real, imagp: &imaginary)
    let length = vDSP_Length(floor(log2(Float(input_windowed.count))))
    let radix = FFTRadix(kFFTRadix2)
    let weights = vDSP_create_fftsetup(length, radix)
    withUnsafeMutablePointer(to: &splitComplex) { splitComplex in
      vDSP_fft_zip(weights!, splitComplex, 1, length, FFTDirection(FFT_FORWARD))
    }
    // zvmags yields power spectrum: |S|^2
    var magnitudes = [Float](repeating: 0.0, count: input_windowed.count)
    withUnsafePointer(to: &splitComplex) { splitComplex in
        magnitudes.withUnsafeMutableBufferPointer { magnitudes in
        vDSP_zvmags(splitComplex, 1, magnitudes.baseAddress!, 1, vDSP_Length(input_windowed.count))
        }
    }
    vDSP_destroy_fftsetup(weights)
    //magnitudes = magnitudes.map{ $0 / Float(input.count) }
    return magnitudes
}

func symmetricZeroPad(signal: inout [Float], length: Int) {
    let originalLength = signal.count
    let zerosToAdd = length - originalLength
    if zerosToAdd <= 0 {
      // No padding needed
      return
    }
    // Pad the signal symmetrically
    let leftZerosCount = zerosToAdd / 2
    let rightZerosCount = zerosToAdd - leftZerosCount
    signal = [Float](repeating: 0.0, count: leftZerosCount) + signal + [Float](repeating: 0.0, count: rightZerosCount)
}

差异原因分析 #

1. 频谱泄漏问题 #

是的，两者的频谱泄漏情况完全不同：

• Python代码未使用任何窗口函数，如果输入的单频正弦波正好是FFT频率点的整数倍（整周期截断），输出会是严格的delta峰，无泄漏。
• Swift代码强制添加了汉明窗，即使是整周期正弦，窗口函数也会让频谱峰展宽，不再是delta函数，这是测试结果差异的核心原因之一。

2. 基函数是否不同 #

两者使用的基函数完全一致：
NumPy的np.fft.fft和Apple vDSP库的vDSP_fft_zip都是基于Radix-2的FFT实现，数学本质相同，都是标准DFT的复指数基函数，不存在基函数差异。

3. 数据传入与预处理差异 #

数据传入本身没有错误，但预处理和后处理步骤存在多处关键差异：

• 零填充：Swift将信号对称零填充到2048点，而Python使用原始信号长度做FFT，导致频率分辨率不同。
• 功率谱计算：
  • Python取FFT结果的前半部分（实信号FFT具有共轭对称性，后半部分冗余），计算(1/N)*|FFT|²后，再归一化除以最大值。
  • Swift保留了完整的FFT长度结果，仅用vDSP_zvmags计算|FFT|²，未做1/N的缩放，也没有归一化。
• FFT长度匹配：Swift中length参数取的是log2(填充后长度)的整数，这部分逻辑正确，但和Python的原始长度FFT不匹配。

修正建议 #

要让Swift输出和Python对齐，需调整以下几点：

1. 移除汉明窗口（或给Python代码添加相同窗口，保持一致）；
2. 取消对称零填充，使用原始信号长度进行FFT；
3. 功率谱计算时，对结果做1/N的缩放，并仅保留前半部分数据；
4. 添加归一化步骤，将结果除以最大值。

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