如何在C#/C++中实现音频匹配所需的scipy.signal.correlate等价功能？

我有一段Python代码，用于在大型音频文件中查找匹配的小型音频片段，核心依赖scipy.signal.correlate函数实现快速互相关。目前已能将其余代码转译为C#，但找不到该函数的等价实现，以下是原Python代码：

import datetime
from scipy import signal
import numpy as np
import librosa

def seconds2str(seconds):
    return str(datetime.timedelta(seconds=seconds))

def find_templ(src_path, templ_path):

    samples_src, sample_rate = librosa.load(src_path, sr=None)
    samples_templ, sample_rate_templ = librosa.load(templ_path, sr=sample_rate)
    src_duration = librosa.get_duration(y=samples_src, sr=sample_rate)
    templ_duration = librosa.get_duration(y=samples_templ, sr=sample_rate_templ)

    if templ_duration > src_duration:
        print("warning: You are looking for a clip within an audio clip shorter than it!!!")

    print(f"src   Duration: {seconds2str(src_duration)}")
    print(f"templ Duration: {seconds2str(templ_duration)}")

    correlate = signal.correlate(samples_src[:sample_rate * int(src_duration)], samples_templ, mode='valid', method='fft')
    start = np.round(np.argmax(correlate) / sample_rate, 2)
    end = start + templ_duration
    print(f"The moment the search stops in seconds: {seconds2str(src_duration)}")
    print(f"Start: {seconds2str(start)}")
    print(f"End:   {seconds2str(end)}")


full_audio = "Track08.mp3"
small_audio = "adv.wav"
find_templ(full_audio, small_audio)

C#实现方案 #

核心思路 #

Python中signal.correlate的method='fft'是通过**快速傅里叶变换(FFT)**实现互相关，避免直接时域运算的高复杂度。C#中可以借助第三方库完成音频读取和FFT运算：

1. 音频读取：用NAudio库加载音频文件，提取PCM样本数据和采样率，替代librosa的功能。
2. 互相关计算：用MathNet.Numerics库的FFT工具实现快速互相关，或者手动实现FFT-based互相关逻辑。

示例代码 #

using System;
using NAudio.Wave;
using MathNet.Numerics;
using MathNet.Numerics.IntegralTransforms;
using System.Numerics;
using System.Linq;

public static class AudioMatcher
{
    private static string SecondsToStr(double seconds)
    {
        return TimeSpan.FromSeconds(seconds).ToString();
    }

    public static void FindTemplate(string srcPath, string templPath)
    {
        // 读取源音频
        float[] srcSamples = LoadAudioSamples(srcPath, out int sampleRate);
        // 读取模板音频，强制匹配源采样率
        float[] templSamples = LoadAudioSamples(templPath, out _, sampleRate);

        double srcDuration = srcSamples.Length / (double)sampleRate;
        double templDuration = templSamples.Length / (double)sampleRate;

        if (templDuration > srcDuration)
        {
            Console.WriteLine("warning: You are looking for a clip within an audio clip shorter than it!!!");
        }

        Console.WriteLine($"src   Duration: {SecondsToStr(srcDuration)}");
        Console.WriteLine($"templ Duration: {SecondsToStr(templDuration)}");

        // 实现FFT-based互相关（对应scipy的method='fft'）
        int fftSize = NextPowerOfTwo(srcSamples.Length + templSamples.Length - 1);
        Complex[] srcFft = new Complex[fftSize];
        Complex[] templFft = new Complex[fftSize];

        // 填充数据到FFT数组
        for (int i = 0; i < srcSamples.Length; i++)
            srcFft[i] = new Complex(srcSamples[i], 0);
        for (int i = 0; i < templSamples.Length; i++)
            templFft[i] = new Complex(templSamples[i], 0);

        // 执行FFT
        Fourier.Forward(srcFft, FourierOptions.Matlab);
        Fourier.Forward(templFft, FourierOptions.Matlab);

        // 计算共轭相乘（互相关的频域等价操作）
        Complex[] crossCorrFft = new Complex[fftSize];
        for (int i = 0; i < fftSize; i++)
            crossCorrFft[i] = srcFft[i] * Complex.Conjugate(templFft[i]);

        // 逆FFT得到时域互相关结果
        Fourier.Inverse(crossCorrFft, FourierOptions.Matlab);

        // 提取valid模式的结果：长度为srcLength - templLength + 1
        int validStart = templSamples.Length - 1;
        int validLength = srcSamples.Length - templSamples.Length + 1;
        double[] corrValues = new double[validLength];
        for (int i = 0; i < validLength; i++)
            corrValues[i] = crossCorrFft[validStart + i].Real;

        // 找到最大值索引
        int maxIndex = Array.IndexOf(corrValues, corrValues.Max());
        double start = Math.Round(maxIndex / (double)sampleRate, 2);
        double end = start + templDuration;

        Console.WriteLine($"The moment the search stops in seconds: {SecondsToStr(srcDuration)}");
        Console.WriteLine($"Start: {SecondsToStr(start)}");
        Console.WriteLine($"End:   {SecondsToStr(end)}");
    }

    private static float[] LoadAudioSamples(string path, out int sampleRate, int targetSampleRate = -1)
    {
        using (var audioFileReader = new AudioFileReader(path))
        {
            sampleRate = audioFileReader.WaveFormat.SampleRate;
            // 如果指定目标采样率，转换采样率
            if (targetSampleRate != -1 && targetSampleRate != sampleRate)
            {
                using (var resampler = new MediaFoundationResampler(audioFileReader, new WaveFormat(targetSampleRate, audioFileReader.WaveFormat.BitsPerSample, audioFileReader.WaveFormat.Channels)))
                {
                    return ReadAllSamples(resampler);
                }
            }
            else
            {
                return ReadAllSamples(audioFileReader);
            }
        }
    }

    private static float[] ReadAllSamples(IWaveProvider reader)
    {
        var samples = new System.Collections.Generic.List<float>();
        var buffer = new float[4096];
        int bytesRead;
        while ((bytesRead = reader.Read(buffer, 0, buffer.Length)) > 0)
        {
            samples.AddRange(buffer.Take(bytesRead));
        }
        return samples.ToArray();
    }

    private static int NextPowerOfTwo(int x)
    {
        x--;
        x |= x >> 1;
        x |= x >> 2;
        x |= x >> 4;
        x |= x >> 8;
        x |= x >> 16;
        return x + 1;
    }

    // 测试调用
    public static void Main()
    {
        string fullAudio = "Track08.mp3";
        string smallAudio = "adv.wav";
        FindTemplate(fullAudio, smallAudio);
    }
}

注意事项 #

• 需要通过NuGet安装NAudio和MathNet.Numerics包。
• 代码中实现的是单声道音频处理，如果是多声道，需要先转成单声道（比如取均值）。

C++实现方案 #

核心思路 #

C++中同样通过FFT实现快速互相关，依赖以下库：

1. 音频读取：用libsndfile库读取各种格式的音频文件，提取PCM样本。
2. FFT运算：用FFTW库实现高效的傅里叶变换。

示例代码 #

#include <iostream>
#include <vector>
#include <cmath>
#include <sndfile.h>
#include <fftw3.h>
#include <algorithm>
#include <cstdio>

std::string SecondsToStr(double seconds)
{
    int hours = static_cast<int>(seconds) / 3600;
    int minutes = (static_cast<int>(seconds) % 3600) / 60;
    double secs = seconds - hours * 3600 - minutes * 60;
    char buffer[64];
    snprintf(buffer, sizeof(buffer), "%02d:%02d:%06.3f", hours, minutes, secs);
    return std::string(buffer);
}

std::vector<float> LoadAudioSamples(const std::string& path, int& sampleRate, int targetSampleRate = -1)
{
    SF_INFO sfInfo;
    SNDFILE* sndFile = sf_open(path.c_str(), SFM_READ, &sfInfo);
    if (!sndFile)
    {
        std::cerr << "Failed to open audio file: " << path << std::endl;
        exit(1);
    }

    sampleRate = sfInfo.samplerate;
    // 读取所有样本，转单声道
    std::vector<float> samples(sfInfo.frames * sfInfo.channels);
    sf_read_float(sndFile, samples.data(), samples.size());
    sf_close(sndFile);

    std::vector<float> monoSamples;
    for (size_t i = 0; i < samples.size(); i += sfInfo.channels)
    {
        monoSamples.push_back(samples[i]);
    }

    // 采样率转换（此处省略，可使用libsamplerate实现）
    if (targetSampleRate != -1 && targetSampleRate != sampleRate)
    {
        std::cerr << "Sample rate conversion not implemented in this example" << std::endl;
        exit(1);
    }

    return monoSamples;
}

void FindTemplate(const std::string& srcPath, const std::string& templPath)
{
    int sampleRate;
    std::vector<float> srcSamples = LoadAudioSamples(srcPath, sampleRate);
    int templSampleRate;
    std::vector<float> templSamples = LoadAudioSamples(templPath, templSampleRate, sampleRate);

    double srcDuration = srcSamples.size() / static_cast<double>(sampleRate);
    double templDuration = templSamples.size() / static_cast<double>(sampleRate);

    if (templDuration > srcDuration)
    {
        std::cout << "warning: You are looking for a clip within an audio clip shorter than it!!!" << std::endl;
    }

    std::cout << "src   Duration: " << SecondsToStr(srcDuration) << std::endl;
    std::cout << "templ Duration: " << SecondsToStr(templDuration) << std::endl;

    // 计算FFT尺寸为2的幂
    int fftSize = 1;
    while (fftSize < srcSamples.size() + templSamples.size() - 1)
        fftSize <<= 1;

    // 分配FFTW内存
    fftw_complex* srcFft = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fftSize);
    fftw_complex* templFft = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fftSize);
    fftw_complex* crossCorrFft = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * fftSize);

    // 初始化数组
    memset(srcFft, 0, sizeof(fftw_complex) * fftSize);
    memset(templFft, 0, sizeof(fftw_complex) * fftSize);
    for (size_t i = 0; i < srcSamples.size(); i++)
        srcFft[i][0] = srcSamples[i];
    for (size_t i = 0; i < templSamples.size(); i++)
        templFft[i][0] = templSamples[i];

    // 创建FFT计划
    fftw_plan srcPlan = fftw_plan_dft_1d(fftSize, srcFft, srcFft, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan templPlan = fftw_plan_dft_1d(fftSize, templFft, templFft, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan invPlan = fftw_plan_dft_1d(fftSize, crossCorrFft, crossCorrFft, FFTW_BACKWARD, FFTW_ESTIMATE);

    // 执行FFT
    fftw_execute(srcPlan);
    fftw_execute(templPlan);

    // 共轭相乘
    for (int i = 0; i < fftSize; i++)
    {
        crossCorrFft[i][0] = srcFft[i][0] * templFft[i][0] + srcFft[i][1] * templFft[i][1];
        crossCorrFft[i][1] = srcFft[i][1] * templFft[i][0] - srcFft[i][0] * templFft[i][1];
    }

    // 逆FFT
    fftw_execute(invPlan);

    // 提取valid模式结果，归一化
    int validStart = templSamples.size() - 1;
    int validLength = srcSamples.size() - templSamples.size() + 1;
    std::vector<double> corrValues(validLength);
    for (int i = 0; i < validLength; i++)
    {
        corrValues[i] = crossCorrFft[validStart + i][0] / fftSize;
    }

    // 找最大值索引
    auto maxIt = std::max_element(corrValues.begin(), corrValues.end());
    int maxIndex = std::distance(corrValues.begin(), maxIt);
    double start = std::round(maxIndex / static_cast<double>(sampleRate) * 100) / 100;
    double end = start + templDuration;

    std::cout << "The moment the search stops in seconds: " << SecondsToStr(srcDuration) << std::endl;
    std::cout << "Start: " << SecondsToStr(start) << std::endl;
    std::cout << "End:   " << SecondsToStr(end) << std::endl;

    // 释放资源
    fftw_destroy_plan(srcPlan);
    fftw_destroy_plan(templPlan);
    fftw_destroy_plan(invPlan);
    fftw_free(srcFft);
    fftw_free(templFft);
    fftw_free(crossCorrFft);
}

int main()
{
    std::string fullAudio = "Track08.mp3";
    std::string smallAudio = "adv.wav";
    FindTemplate(fullAudio, smallAudio);
    return 0;
}

注意事项 #

• 需要安装libsndfile和FFTW库，编译时链接对应的库文件。
• 示例中省略了采样率转换逻辑，可通过libsamplerate库补充实现。
• 多声道音频需先转换为单声道再处理。

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