Ускорение и NumPy дают разные результаты для БПФ - PullRequest
1 голос
/ 12 марта 2019

Я работаю над конвейером разработки функций для CoreML, и мне нужно выполнить БПФ для моих данных. Проблема заключается в том, что результаты инфраструктуры ускорения и результаты БПФ NumPy различны.

Swift:

public func fft(_ input: [Double]) -> [Double] {
    var real = [Double](input)
    var imaginary = [Double](repeating: 0.0, count: input.count)
    var splitComplex = DSPDoubleSplitComplex(realp: &real, imagp: &imaginary)

    let length = vDSP_Length(floor(log2(Float(input.count))))
    let radix = FFTRadix(kFFTRadix2)
    let weights = vDSP_create_fftsetupD(length, radix)
    vDSP_fft_zipD(weights!, &splitComplex, 1, length, FFTDirection(FFT_FORWARD))

    var magnitudes = [Double](repeating: 0.0, count: input.count)
    vDSP_zvmagsD(&splitComplex, 1, &magnitudes, 1, vDSP_Length(input.count))

    var normalizedMagnitudes = [Double](repeating: 0.0, count: input.count)
    vDSP_vsmulD(sqrt(magnitudes), 1, [2.0 / Double(input.count)], &normalizedMagnitudes, 1, vDSP_Length(input.count))

    vDSP_destroy_fftsetupD(weights)

    return normalizedMagnitudes
}

Python:

def fft(series: pd.Series):
    f = np.fft.fft(series)
    fa = np.abs(f)
    return pd.Series(fa)

Я использую одинаковые 100 значений для каждого метода.

Полагаю, это как-то связано с нормализующей частью, но я даже не уверен, что оба массива содержат одно и то же:

  • индекс 0: нулевой частотный член
  • индекс 1-50: положительные величины
  • индекс 50-99: отрицательные величины

Меня интересуют только положительные величины.

Редактировать:

Вот график NumPy:

NumPy

А вот и Ускоренный сюжет:

Accelerate

Я надеюсь, что кто-то может помочь:)

Ответы [ 2 ]

2 голосов
/ 12 марта 2019

Есть две проблемы:

  • Как указано E.Coms , реализация с использованием FFT инфраструктуры Accelerate включает в себя шаг нормализации, который берет квадратный корень из величины и умножается на скаляр 2/N. Реализация, использующая NumPy, отсутствует.

  • БПФ NumPy поддерживает входы произвольной длины, и результирующие частотные интервалы соответствуют ожидаемым (нулевая частота в индексе 0, положительные частоты в индексе 1-50 и отрицательные частоты в индексе 51-99). С другой стороны, FFT в платформе Accelerate должна иметь длину, равную степени 2. Соответственно, этот пример кода вычисляет FFT из ваших первых 64 входных значений. Это помещает нулевую частоту в индекс 0, положительные частоты в индекс 1-32 и отрицательные частоты в индекс 33-63. Остальные индексы (64-99) являются просто оставшимися нетронутыми входами.

0 голосов
/ 12 марта 2019

Если ускорение возвращает квадратную величину, результат будет таким же, как python.

  public func fft(_ input: [Double]) -> [Double] {
   ....
   return  magnitudes.map{sqrt($0)}
   }

Вы можете выяснить, какой из них не сейчас.Если вы хотите использовать ускорение, можете использовать следующее:

    public func fft(_ input: [Double]) -> [Double] {
   ....

    var normalizedMagnitudes = [Double](repeating: 0.0, count: input.count)
    var count : Int32 = Int32(input.count)
    vvsqrt(   &normalizedMagnitudes, &magnitudes, &count )



     vDSP_destroy_fftsetupD(weights)

     return normalizedMagnitudes}
...