Анимируем свёртку функций

Проведём пошаговую операцию свёртки двух сигналов, чтобы детально разобраться в смысле данной математической операции.

Свёртка это математическая операция, применимая к и , которая порождает третью функцию и определяется формулой:

Рассмотрим реализацию свёртки двух сигналов на примере. Пусть сигналы и задаются в следующем виде:

Их графики, соответственно будут представлять собой ломаные, как показано ниже

Графики сигналов f(t) и g(t)

Пример свёртки двух функций пошагово

Итак, начнём с того, что в интегральном выражении у функции по переменной, которой производится интегрирование, стоит знак "-", а значит график нужно отразить по горизонтальной оси относительно оси ординат, как показано на рисунке ниже:

Отражённый сигнал g(t)

Выражение означает, что движение будет происходить оси . Причём, если , то график функции будет двигаться влево, при движение будет происходить вправо.

Далее необходимо вычислить произведение для каждого значения и посчитать площадь получившейся в результате произведения двух графиков фигуры.

На рисунке ниже показано, какие случаи пересечений нам необходимо просчитать по мере движения графика по оси .

Графичекая пошаговая интерпретация свёртки сигналов f(t) и g(t)

Остановимся на каждом из 6-ти случаев и произведём расчёты площадей получившихся фигур.

При графики сигналов не пересекаются, а значит и площадь фигуры их пересечения нулевая, то есть, результирующий интеграл примет вид:
На интервале сигнал представляет собой константу: , а график сигнала - есть прямая . Исходя из вышесказанного, результатом сворачивания двух функций будет функция вида:
На интервале результирующий интеграл можно разбить на 2 интеграла: первый вида, описанного в пункте b). Второй интеграл будет представлять собой площадь фигуры, имеющей стороны, описываемые прямыми и . В результате, функция примет вид:
На интервале результатом свёртки будет функция :
Предпоследним, дающим непустое пересечение графиков, будет интервал , который породит функцию :
Результирующей функцией на интервале будет второе слагаемое из предыдущего случая со своими пределами интегрирования:

Таким образом, результатом свёртки игналов и будет функция , задаваемая совокупностью уравнений:

На рисунке ниже приведена функция являющаяся результатом сворачивания двух сигналов на каждом из смысловых промежутков.

Функция y(t) - результат свёртки сигналов f(t) и g(t)

Реализация свёртки на Scilab

Для реализации в численных методах, рассматривается дискретный случай свертки двух функций и

где - количество элементов функций и соответственно.

Прежде всего, зададим функции и с помощью условных операторов в Scilab:

function ft = f(t)
    ft = zeros(t);
    
    for j = 1:length(t)   
        if ( t(j) >= 0 && t(j) < 1 ) then
            ft(j) = 2;
        elseif ( t(j) >= 1 && t(j) < 2 )  then
            ft(j) = -2*t(j) + 4;  
        end;
    end;
endfunction

Задание сигнала f(t)

function gt = g(t)  
    gt = zeros(t);
    
    for j = 1:length(t)   
        if ( t(j) == 0 ) then
            gt(j) = 1;
        elseif ( t(j) >= 0 && t(j) < 3 )  then
            gt(j) = -1/3*t(j) + 1;  
        end;
    end;
endfunction

Задание сигнала g(t)

Далее на пишем функцию, реализующую свёртку, согласно формулы (3)

function cnv = myconv(f_, g_)
    Lf = length(f_);
    Lg = length(g_);
    
    L = Lf + Lg - 1;
    cnv = zeros(1,L);
    
    for k = 1:L 
        jmin = max(1, k + 1 - Lg);
        jmax = min(k, Lf);
      
        for j = jmin:jmax      
            cnv(k) = cnv(k) + f_(j) * g_(k + 1 - j);
        end        
    end
endfunction

Задание сигнала g(t)

И посмотрим, сойдётся ли результат, полученный с помощью нашей функции с результатом, полученным при использовании встроенной функции

d = 1e-1;
Tmin = -1; Tmax = 4;
t = Tmin:d:Tmax;
  
plot( conv(f(t), g(t)), 'r' );
plot( myconv(f(t), g(t)), 'b--' );
xgrid();xtitle("Свёртка двух сигналов","t", "y(t)");
legend("Встроенная свёртка", "Своя свёртка");
gca().children.children.thickness  = 2;

Вывод графиков двух функций с настройкой толщины в Scilab

Глядя на изображение ниже, можно с уверенностью сказать, что результаты совпадают :)

Результат наложения встроенной и самописной функции свёртки

Вот только масштаб и отсчёты не похожи на наши задынные промежутки. Поэтому, проведём масштабирование результатов.

Корректировка результатов свёртки

Для корректного отображения результирующего графика нам потребуется несколько модифицировать функцию . Прежде всего, заведём переменную, отвечающую за сдвиг результата по горизонтальной оси:

N = (2*abs(Tmin))/d;

Переменная сдвига

И изменим функцию следующим образом:

function cnv = myconv(f_, g_, shift)
    Lf = length(f_);
    Lg = length(g_);
    
    L = Lf + Lg - 1;
    cnv = zeros(2,L);
    
    for k = 1:L 
        jmin = max(1, k + 1 - Lg);
        jmax = min(k, Lf);
        cnv(1,k) = k - shift;
        for j = jmin:jmax      
            cnv(2,k) = cnv(2,k) + f_(j) * g_(k+1-j);
        end        
    end
endfunction

Обновлённая функция свёртки

Посмотрим, как теперь будет выглядеть результат:

myconv_scalled = d * myconv(f(t), g(t), N);
plot(myconv_scalled(1,:), myconv_scalled(2,:), 'b--' );
xgrid();xtitle("Свёртка двух сигналов","t", "y(t)");
legend("Своя свёртка");
gca().children.children.thickness  = 2;

Рисум график масштабированной свёртки

Результат модификации свёртки

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

Реализация анимации свёртки двух сигналов

Прежде всего, выведем графики сигналов, которые сворачиваем:

shift = 2*abs(Tmin);
plot(t, f(t), "r" ); 
plot(-t - shift, g(t), "g" );
xgrid();xtitle("Свёртка двух сигналов","t", "f(t), g(t), y(t)");
legend("f(t)", "g(-t)");

conv_axes = gca();
conv_axes.data_bounds=[-5 -0.1; 8 max(myconv_scalled(2,:))];
conv_axes.children.children.thickness  = 2;

Рисуем графики функций f(t) и отражённой и сдвинутой g(t)

Чтобы графики выглядели поприличнее, нужно уменьшить шаг дискретизации, например до и настроить их толщину с помощью

Графики функций f(t) и g(t)

Осталось лишь реализовать анимационный цикл.

В статье дано описание основ анимации в Scilab.

Здесь на каждом шагу цикла мы рисуем i-ю точку свёртки, а зелёный график рисуем и стираем, что имитирует его движение.

i = 1; t_ = Tmin; g_inv = g(t);
while ( (i <= length(myconv_scalled(2,:))) && (t_ < 2*Tmax)  )
    drawlater();    
        sca(conv_axes);           
                      
        plot(myconv_scalled(1,i), myconv_scalled(2,i), '.-b');
        plot(-t - shift + i*d, g_inv, 'g');
        conv_axes.children.children.thickness  = 2;
        
        realtime(i);
        delete(conv_axes.children.children(3));
    drawnow();    
    
   i = i+1;
   t_ = i*d;
end