scilab

//BcIb4
//порядок 4
//фильтр баттерворта - каскадная схема из прямых форм
//фильтр Чебышева II рода, прямая форма II (каноническая)
//_______________________________________________________
//Модельный сигнал
N=1024;
dT=0.01;
x=0:dT:N*dT;
sig=0.7*sin(2*%pi*2*x)+0.35*cos(2*%pi*4*x);
y=sig+grand(1, (N+1), "nor", 0, 0.2);
clf(0); scf(0);
subplot(211);
plot2d(x,y);
xtitle("График модельного сигнала","Время t, с", "Амплитуда, ед.");
//Спектр модельного сигнала
spectr = abs(fft(y));
dW = 1/(dT*N); //Шаг по частота
W=0:dW:1/dT; //Диапазон частот
_W_ = length(W);
_spectr_=length(spectr);
subplot(212);
plot2d(W(1:_W_/2),spectr(1:_spectr_/2)); //Делим для вывода спектра до частоты Найквиста
xtitle("Спектр сигнала","Частота f, Гц","Амплитуда, ед.");

//Аналоговые фильтры
hsButt=analpf(4, 'butt', [0 0], 0.3*2/dT); //строим вид передаточной функции для фильров
hsCheb2=analpf(4,'cheb2',[0.1 0.1],0.3*2/dT);
ls_Butt=syslin('c', hsButt); //задаём линейную систему
ls_Cheb2=syslin('c', hsCheb2);
scf(1);clf(1);
bode(ls_Butt); //Диаграмма Боде для полученной линейной системы
xtitle("АЧХ и ФЧХ фильтра Баттерворта");
scf(2);clf(2);
bode(ls_Cheb2);
xtitle("АЧХ и ФЧХ фильтра Чебышёва II рода");
//Обработка спектра сигнала фильтрами
FButt=repfreq(ls_Butt,W);
FCheb2=repfreq(ls_Cheb2,W);
SpectrButt=fft(y).*FButt;//Обработка спектра
SpectrCheb2=fft(y).*FCheb2;
scf(3);clf(3);
subplot(311);
plot2d(W(1:_W_/2),spectr(1:_spectr_/2));
xtitle("Спектр сигнала до фильтрации","Частота f, Гц", "Aмплитуда, ед.");
subplot(312);
plot2d(W(1:_W_/2),abs(SpectrButt(1:_spectr_/2)));
xtitle("Спектр сигнала после фильтра Баттерворта","Частота f, Гц", "Aмплитуда, ед.");
subplot(313);
plot2d(W(1:_W_/2),abs(SpectrCheb2(1:_W_/2)));
xtitle("Спектр сигнала после фильтра Чебышёва II рода","Частота f, Гц", "Aмплитуда, ед.");
//Обратноe прeобразованиe обработанного спeктра
reverseButt=real(ifft(SpectrButt)); //Восстановление из спектра
reverseCheb2=real(ifft(SpectrCheb2));
reverseButt=reverseButt/max(reverseButt); //Нормировка по максимальному значению
reverseCheb2=reverseCheb2/max(reverseCheb2);
clf(4);scf(4);
subplot(311);
plot2d(x, (y/max(abs(y))));
xtitle("Сигнал до фильтрации","Время t, с", "Амплитуда, ед.");
subplot(312);
plot2d(x, reverseButt);
xtitle("Сигнал после фильтра Баттерворта","Время t, с", "Амплитуда, ед.");
subplot(313);
plot2d(x, reverseCheb2);
xtitle("Сигнал после фильтра Чебышёва II рода","Время t, с", "Амплитуда, ед.");
//Переходим к цифровым фильтрам
D_hsButt=iir(4,'lp','butt',[0.3 0], [0 0]);
//                                      2            3            4
//   0.1671793 + 0.6687171z + 1.0030756z + 0.6687171z + 0.1671793z
//   --------------------------------------------------------------
//                                          2            3   4
//       0.0301189 + 0.1826757z + 0.6799785z + 0.7820952z + z
D_hsCheb2=iir(4,'lp','cheb2',[0.3 0], [0.1 0.1]);
//                                      2            3            4
//   0.2727433 + 0.6738200z + 0.8994179z + 0.6738200z + 0.2727433z
//   --------------------------------------------------------------
//                                          2            3   4
//       0.0780474 + 0.2898983z + 0.8045473z + 0.6200515z + z
[F_Butt,D_butt]=repfreq(D_hsButt,0,0.5);
clf(5);scf(5);
subplot(211);
plot2d(F_Butt*100,abs(D_butt));
xtitle("АЧХ дискретного фильтра Баттерворта","Частота f, Гц", "Амплитуда");
//https://help.scilab.org/doc/5.5.2/en_US/dbphi.html
[db,phi]=dbphi(D_butt);
subplot(212);
plot2d(F_Butt*100, phi);
xtitle("ФЧХ дискретного фильтра Баттерворта","Частота f, Гц", "Фаза");
[F_Cheb2,D_Cheb2]=repfreq(D_hsCheb2,0,0.5);
clf(6);scf(6);
subplot(211);
plot2d(F_Cheb2*100,abs(D_Cheb2));
xtitle("АЧХ дискретного фильтра Чебышёва II рода","Частота f, Гц", "Амплитуда");
[db,phi]=dbphi(D_Cheb2);
subplot(212);
plot2d(F_Cheb2*100,phi);
xtitle("ФЧХ дискретного фильтра Чебышёва II рода","Частота f, Гц)", "Фаза");


//каскадная схема из прямых форм для баттерфорда
//вторая строка коэффициентов с обратным знаком
//в соответсвии со схемой - учитываем знаки и берем обратные
K=factors(D_hsButt(3));
K_1=coeff(c(1));
K_2=coeff(c(2));
D_ButtGr=y;
D_ButtGr2=y;
k1=1;
k2=1;
for i=3:N
    D_ButtGr(i)=k1*y(i)+2*k1*y(i-1)+k1*y(i-2)-K_1(2)*D_ButtGr(i-1)-K_1(1)*D_ButtGr(i-2);
    D_ButtGr2(i)=k2*D_ButtGr(i)+2*k2*D_ButtGr(i-1)+k2*D_ButtGr(i-2)-K_2(2)*D_ButtGr(i-1)-K_2(1)*D_ButtGr2(i-2);
end

//прямая схема (прямая форма II) каноническая для фильтра Чебышева II рода
//определение числителя и знаменателя
a = coeff(D_hsCheb2(3));
b = coeff(D_hsCheb2(2));
//
DCHEB = zeros(1,length(y)+5);
v = zeros(1,length(y)+5);
//реализуем схему
for i=6:N
    v(i) = y(i) - (a(4)*v(i-1) + a(3)*v(i-2) + a(2)*v(i-3) + a(1)*v(i-4));
    DCHEB(i) = b(5)*v(i-1) + b(4)*v(i-2) + b(3)*v(i-3) +b(2)*v(i-4) + b(1)*v(i-5);
end

G_DCHEB = DCHEB; //Временная переменная для графика, где-то что-то переопределяю, поэтому не спроилось, костыль :(

D_ButtGr2=D_ButtGr2/max(D_ButtGr2);
MDCHEB=DCHEB/max(DCHEB);
//Спектр сигнала
SpecD_ButtGr=fft(D_ButtGr2);
SpMDCHEB = fft(MDCHEB);
clf(7);scf(7);
subplot(311);
plot2d(W(1:_W_/2),spectr(1:_spectr_/2));
xtitle("Спектр сигнала до фильтрации","Частота f, Гц", "Амплитуда, ед.");
subplot(312);
plot2d(W(1:_W_/2),abs(SpecD_ButtGr(1:_W_/2)));
xtitle("Спектр сигнала после фильтра Баттерворта каскадной формы","Частота f, Гц", "Амплитуда, ед.");
subplot(313);
plot2d(W(1:_W_/2),abs(SpMDCHEB(1:_W_/2)));
xtitle("Спектр сигнала после фильтра Чебышева II рода","Частота f, Гц", "Амплитуда, ед.");

//Сигнал после фильтрации
clf(8);scf(8);
subplot(311);
plot2d(x,(y/max(abs(y))));
xtitle("Исходный сигнал","Время t, с", "Амплитуда, ед.");
subplot(312);
plot2d(x(1:length(D_ButtGr2)),D_ButtGr2);
xtitle("Сигнал после обработки цифровым фильтром Баттерворта каскадной формы","Время t, с", "Амплитуда, ед.");
subplot(313);
plot2d(x,G_DCHEB(6:length(G_DCHEB)));
xtitle("Сигнал после обработки цифровым фильтром Чебышева","Время t, с", "Амплитуда, ед.");