Есть такой способ: Можно подобрать цифровой фильтр с подобной для низких частот частотной характеристикой и програмно
реализовать его. Пусть у
вас есть цифровая система которая работает с периодом Δt. В случае
низкочастотных сигналов Δt << k (k в вашей формуле,
обычно для параметра используют символ Т), достаточно заменить s
=(2/Δt)(1-z-1)/(1+z-1)
И получить системную функцию фильтра: W(s)~ H(z) =
1/(1+ks)=1/(1+k(2/Δt)(1-z-1)/(1+z-1)) = A(1+z-1)/(1-B z-1)
где A = 1/(1+2k/Δt) , B = (2k/Δt-1)/(1+2k/Δt)
Дальше выражаете yn и готово.
Фильтр будет смещать в сторону меньших частот нули и полюсы сигналов
ωц = (2/Δt)arctg(ωаΔt/2). Если это критично можно сначала ввести предискажение на нужных частотах, а потом преобразовывать
Какое-то видео как это делать