如何用插值将double[]数组降采样至任意长度？

嘿，针对你这个FFT结果降采样的需求，我刚好有相关实践经验，咱们一步步拆解解决：

一、先明确FFT数据的核心特性，避免踩坑 #

FFT数组里的每个double值对应一个频率bin——它代表的是一段连续频率区间的能量积分，而非单点频率。直接对bin值做简单平均会丢失峰值细节，而插值法刚好能满足你"任意输出长度、保留极端电平、无舍入问题"的需求。

二、基于Math.NET Numerics的具体实现 #

1. 先构建频率轴（关键前提） #

你需要先明确原FFT数组对应的频率范围，假设原数组长度为4096，原始信号采样率为sampleRate（单位Hz），可以这样生成频率轴：

// 双边FFT的频率轴（范围0到sampleRate）
double[] originalFrequencies = Enumerable.Range(0, 4096)
    .Select(i => (double)i * sampleRate / 4096)
    .ToArray();

// 如果是单边FFT（仅正频率，长度2048），则频率范围到sampleRate/2
// double[] originalFrequencies = Enumerable.Range(0, 2048)
//     .Select(i => (double)i * sampleRate / 4096)
//     .ToArray();

2. 用CubicSpline实现平滑插值（满足多源值组合需求） #

Math.NET的CubicSpline刚好符合你的要求：它基于相邻4个源点计算加权组合值，既能保证平滑性，又不会丢失峰值细节。核心代码如下：

using MathNet.Numerics.Interpolation;

public double[] ResampleFftSpectrum(double[] originalAmplitudes, double sampleRate, int targetLength)
{
    int originalLength = originalAmplitudes.Length;
    
    // 生成原频率轴
    double[] originalFreqs = Enumerable.Range(0, originalLength)
        .Select(i => (double)i * sampleRate / originalLength)
        .ToArray();
    
    // 生成目标频率轴：和原频率范围一致，均匀分布为targetLength个点
    double[] targetFreqs = Enumerable.Range(0, targetLength)
        .Select(i => (double)i * sampleRate / targetLength)
        .ToArray();
    
    // 创建立方样条插值器（要求输入的频率轴已排序，这里天然满足）
    var spline = CubicSpline.InterpolateSorted(originalFreqs, originalAmplitudes);
    
    // 在目标频率轴上计算插值结果
    return targetFreqs.Select(f => spline.Interpolate(f)).ToArray();
}

3. 优化：强制保留极端峰值电平 #

如果你的FFT结果存在尖锐的单频率峰值（比如某个bin幅值远高于周围），可以先提取峰值点，再基于这些点做插值，确保峰值不会被平滑掉：

public double[] ResampleFftWithPeakPreservation(double[] originalAmplitudes, double sampleRate, int targetLength)
{
    int originalLength = originalAmplitudes.Length;
    double[] originalFreqs = Enumerable.Range(0, originalLength)
        .Select(i => (double)i * sampleRate / originalLength)
        .ToArray();
    
    // 提取峰值点（简单检测：比左右相邻点都高）
    var peakIndices = new List<int> { 0, originalLength - 1 }; // 先加入首尾点
    for (int i = 1; i < originalLength - 1; i++)
    {
        if (originalAmplitudes[i] > originalAmplitudes[i-1] && originalAmplitudes[i] > originalAmplitudes[i+1])
        {
            peakIndices.Add(i);
        }
    }
    peakIndices.Sort(); // 确保顺序正确
    
    // 构建包含峰值的插值点集
    var peakFreqs = peakIndices.Select(i => originalFreqs[i]).ToArray();
    var peakAmplitudes = peakIndices.Select(i => originalAmplitudes[i]).ToArray();
    
    // 基于峰值点做立方样条插值
    var spline = CubicSpline.InterpolateSorted(peakFreqs, peakAmplitudes);
    
    double[] targetFreqs = Enumerable.Range(0, targetLength)
        .Select(i => (double)i * sampleRate / targetLength)
        .ToArray();
    
    return targetFreqs.Select(f => spline.Interpolate(f)).ToArray();
}

三、解决舍入与混叠问题 #

• 舍入问题：插值法基于连续频率轴计算，不管目标长度和原长度是否整除，每个输出值都是精确的加权组合，完全避免了分组平均时的组大小不一致问题。
• 混叠问题：因为我们是在频率域直接重采样，而非时域降采样，只要插值方法合理（比如立方样条），就不会引入混叠——原FFT数据本身已经是频率域的离散表示，插值只是对连续频率曲线的采样。

四、为什么不建议用位图或音频库？ #

• 位图方法：需要将double数组转成像素格式（比如byte），会丢失精度，而且数据转换过程额外消耗性能，完全没必要在数值域能直接处理的情况下绕路。
• 音频库：音频重采样是针对时域信号的，虽然原理类似，但音频库的API并不适配频率域的FFT数据，反而会增加不必要的依赖，不如直接用你已经在使用的Math.NET Numerics。

额外建议：处理复数FFT结果 #

如果你需要保留相位信息，建议分别对FFT的实部和虚部做插值，再计算幅值和相位，比直接插值幅值更准确：

public (double[] resampledAmplitudes, double[] resampledPhases) ResampleComplexFft(
    double[] originalReal, double[] originalImaginary, double sampleRate, int targetLength)
{
    int originalLength = originalReal.Length;
    double[] originalFreqs = Enumerable.Range(0, originalLength)
        .Select(i => (double)i * sampleRate / originalLength)
        .ToArray();
    
    double[] targetFreqs = Enumerable.Range(0, targetLength)
        .Select(i => (double)i * sampleRate / targetLength)
        .ToArray();
    
    // 分别插值实部和虚部
    var splineReal = CubicSpline.InterpolateSorted(originalFreqs, originalReal);
    var splineImag = CubicSpline.InterpolateSorted(originalFreqs, originalImaginary);
    
    double[] resampledReal = targetFreqs.Select(f => splineReal.Interpolate(f)).ToArray();
    double[] resampledImag = targetFreqs.Select(f => splineImag.Interpolate(f)).ToArray();
    
    // 计算幅值和相位
    double[] resampledAmplitudes = resampledReal.Zip(resampledImag, (r, i) => Math.Sqrt(r*r + i*i)).ToArray();
    double[] resampledPhases = resampledReal.Zip(resampledImag, (r, i) => Math.Atan2(i, r)).ToArray();
    
    return (resampledAmplitudes, resampledPhases);
}

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