Реализация фильтра Фарроу на си

[Rem]

Здравствуйте. Я решил написать функцию для использования в программах.
Статья и теория тут - http://www.dsplib.ru/content/farrow/farrow.html
В функцию передают 4 точки и искомую координату. Закомментил часть, где вручную можно вводить значения. Для образца я взял параметры из статьи. Циклом вычисляю 10 значений из среднего диапазона.
Просьба к форумчанам проверить работоспособность.

Код: Выделить всё

#include <stdio.h>
#include <stdlib.h>

float Farrow_Filter(float x1,float y1,float x2,float y2,float x3, float y3, float x4, float y4, float x_current);

int main()
{
    printf("Farrow filter\n");

    float i, x1,y1,x2,y2,x3,y3,x4,y4,x_current, y_current; // 4 точки и искомая координата

    /*    printf("Input 4 points\n");
        printf ( "x1= " ), scanf ("%f",&x1);
        printf ( "y1= " ), scanf ("%f",&y1);
        printf ( "x2= " ), scanf ("%f",&x2);
        printf ( "y2= " ), scanf ("%f",&y2);
        printf ( "x3= " ), scanf ("%f",&x3);
        printf ( "y3= " ), scanf ("%f",&y3);
        printf ( "x4= " ), scanf ("%f",&x4);
        printf ( "y4= " ), scanf ("%f",&y4);
        printf ( "x_current= %f/n" ,x_current);

        */
    x1=0;
    y1=1;
    x2=0.001;
    y2=-1;
    x3=0.002;
    y3=2;
    x4=0.003;
    y4=1;
    x_current=0.0013;

    for (i=0.001; i<0.002; i+=0.0001)
    {
        y_current=Farrow_Filter(x1,y1,x2,y2,x3,y3,x4,y4,i);
        printf("y(%f)=%f \n",i, y_current);
    }
    return 0;
}

float Farrow_Filter(float t0,float y0,float t1,float y1,float t2, float y2, float t3, float y3, float t_current)
{
    //модифицированный фильтр Фарроу http://www.dsplib.ru/content/farrow/farrow.html
    float x,a0,a1,a2,a3,y;

    // нормировка
    x=(t_current-t0)*3/(t3-t0)-2;
    // коэффициенты полинома
    a3=(y3-y0)/6 + (y1-y2)/2;
    a1=(y3-y1)/2 - a3;
    a2=y3-y2-a1-a3;
    a0=y2;
    // вычисляем значение в заданой точке
    y=((x*a3+a2)*x+a1)*x+a0;
    return(y);
}

[Бахурин Сергей]

Выглядит рабочим. Данные выдает адекватные. Точность данных я не проверял потому как функция написана в формате float (не double)

[Rem]

Подскажите, как можно используя эти сплайны решить пару задач?
1. Известно, что в центральном отрезке сплайн пересекает ось. Найти X координату пересечения.
2. Известно, что в центральном отрезке сплайн выгибается в виде холма. Нужно найти X координату вершины.

[Бахурин Сергей]

Ну у вас по сути есть полином


Чтобы найти ноль достатчоно решить это кубическое уравнение. Нахождение максимума это производная полинома должна быть равна нулю. Т.е. надо решить квадратное уравнение.

[Rem]

Сделал поиск вершины. Сама функция:

Код: Выделить всё

float Find_Peak_Farrow(float t0,float y0,float t1,float y1,float t2, float y2, float t3, float y3)
{
    // поиск максимума сплайна
    float a0, a1, a2, a3, y, xmax;
    float A, B, C;
    float d, x1, x2, x_1, x_2, y_1, y_2;

    // коэффициенты полинома
    a3=(y3-y0)/6 + (y1-y2)/2;
    a1=(y3-y1)/2 - a3;
    a2=y3-y2-a1-a3;
    a0=y2;

    printf("\nAx^3 + Bx^2 + Cx + D = 0\n");
    printf("A=%f\n",a3);
    printf("B=%f\n",a2);
    printf("C=%f\n",a1);
    printf("D=%f\n",a0);
    // берём производную, преобразуем в квадратное уравнение
    A=3*a3;
    B=2*a2;
    C=a1;
    printf("\nAx^2 + Bx + C = 0\n");
    printf("A=%f\n",A);
    printf("B=%f\n",B);
    printf("C=%f\n\n",C);

    // решаем квадратное уравнение
    d=B*B-4*A*C;
    printf("d=%f\n",d);
    if (d>0)
    {
        x1=(-B-sqrt(d))/(2*A);
        x2=(-B+sqrt(d))/(2*A);
    }
    printf("x1=%f\n",x1);
    printf("x2=%f\n",x2);

    // преобразуем в масштаб шкалы
    x_1= ((x1+2)*(t3-t0) + 3*t0 )/3.; // точки экстремума
    x_2= ((x2+2)*(t3-t0) + 3*t0 )/3.;
    printf("x_1=%f\n",x_1);
    printf("x_2=%f\n",x_2);

    //считаем значение полинома в точках локальных экстремумах
    y_1 = Farrow_Filter(t0,y0,t1,y1,t2,y2,t3,y3,x_1);
    y_2 = Farrow_Filter(t0,y0,t1,y1,t2,y2,t3,y3,x_2);
    printf("y_1=%f\n",y_1);
    printf("y_2=%f\n",y_2);

    if (y_2 > y_1) x_1 = x_2; // возвращаем бОльший из экстремумов
    return(x_1);
}

Алгоритм:
Вычисляем коэффициенты полинома, берём производную для определения коэффициентов квадратного уравнения, находим две точки экстремумов (где касательная горизонтальна), вычисляем значение функции в каждом, берём максимальное, возвращаем координату при которой есть максимум. Для проверки вычисляем значение функции в этой найденной точке и двух близко расположенных.
Код основной функции может быть такой:

Код: Выделить всё

 int main()
{
    printf("Farrow filter\n");

    float i, x1,y1,x2,y2,x3,y3,x4,y4, y_current, x_max;
    x1=11;
    y1=12;
    x2=12;
    y2=28;
    x3=13;
    y3=31;
    x4=14;
    y4=17;    
    
    x_max= Find_Peak_Farrow(x1,y1,x2,y2,x3,y3,x4,y4); // найти вершину максимума
    printf("x_max=%f \n",x_max);
    y_current=Farrow_Filter(x1,y1,x2,y2,x3,y3,x4,y4,x_max);
    printf("y(%f)=%f \n",x_max, y_current);
    y_current=Farrow_Filter(x1,y1,x2,y2,x3,y3,x4,y4,x_max-0.005);
    printf("y(%f-0.005)=%f \n",x_max, y_current);
    y_current=Farrow_Filter(x1,y1,x2,y2,x3,y3,x4,y4,x_max+0.005);
    printf("y(%f+0.005)=%f \n",x_max, y_current);

    return 0;
}

Для наглядности выводил расчёты на консоль:

Сейчас можно искать вершины волны в звуке и измерять расстояние между ними, вычислять период и частоту.

[ivan219]

Если вы хотите найти длину периода то рекомендую все же искать не через вершины а через 0.