Как вы выполняете бикубическую (или другую нелинейную) интерполяцию передискретизированных аудиоданных?

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

            int newlength = (int)Math.Round(rawdata.Length * lengthMultiplier);
            float[] output = new float[newlength];

            for (int i = 0; i < newlength; i++)
            {
                float realPos = i / lengthMultiplier;
                int iLow = (int)realPos;
                int iHigh = iLow + 1;
                float remainder = realPos - (float)iLow;

                float lowval = 0;
                float highval = 0;
                if ((iLow >= 0) && (iLow < rawdata.Length))
                {
                    lowval = rawdata[iLow];
                }
                if ((iHigh >= 0) && (iHigh < rawdata.Length))
                {
                    highval = rawdata[iHigh];
                }

                output[i] = (highval * remainder) + (lowval * (1 - remainder));
            }

Это работает отлично, но он звучит нормально, только когда я понижаю частоту воспроизведения (т.е. замедляю его). Если я подниму шаг при воспроизведении, этот метод имеет тенденцию производить высокочастотные артефакты, по-видимому, из-за потери информации о образцах.

Я знаю, что бикубические и другие методы интерполяции повторяются с использованием не только двух ближайших значений образца, как в моем примере кода, но я не могу найти хорошие образцы кода (предпочтительно С#), которые я мог бы подключить, чтобы заменить мою линейную метод интерполяции.

Кто-нибудь знает какие-либо хорошие примеры, или кто-нибудь может написать простой метод бикубической интерполяции? Я сделаю это, если придется.:)

Обновление: вот пара реализаций интерполяции С# (благодаря Донни ДеБоеру для первого и ноздредна для второго):

    public static float InterpolateCubic(float x0, float x1, float x2, float x3, float t)
    {
        float a0, a1, a2, a3;
        a0 = x3 - x2 - x0 + x1;
        a1 = x0 - x1 - a0;
        a2 = x2 - x0;
        a3 = x1;
        return (a0 * (t * t * t)) + (a1 * (t * t)) + (a2 * t) + (a3);
    }

    public static float InterpolateHermite4pt3oX(float x0, float x1, float x2, float x3, float t)
    {
        float c0 = x1;
        float c1 = .5F * (x2 - x0);
        float c2 = x0 - (2.5F * x1) + (2 * x2) - (.5F * x3);
        float c3 = (.5F * (x3 - x0)) + (1.5F * (x1 - x2));
        return (((((c3 * t) + c2) * t) + c1) * t) + c0;
    }

В этих функциях x1 представляет собой примерное значение перед точкой, которую вы пытаетесь оценить, а x2 - это значение выборки после вашей точки. x0 слева от x1, а x3 - справа от x2. t идет от 0 до 1 и является расстоянием между оценкой точки и точкой x1.

Метод Эрмита, похоже, работает очень хорошо и, похоже, несколько уменьшает шум. Что еще более важно, кажется, что звук лучше, когда волна ускоряется.

Ответы

double InterpCubic(double x0, double x1, double x2, double x3, double t)
{
   double a0, a1, a2, a3;

   a0 = x3 - x2 - x0 + x1;
   a1 = x0 - x1 - a0;
   a2 = x2 - x0;
   a3 = x1;

   return a0*(t^3) + a1*(t^2) + a2*t + a3;
}