% исходная частота сигнала
f = 90000;
% период частоты
T = 1/f;
% частота выборок 360 кГц
fs  = 360000;
% начало модели
tmin = 0;
% конец модели
tmax = 5*T;
% шаг моделирования
dt = T/4;
% цикл моделирования
t = tmin:dt:tmax;
xc = cos(2*pi*f*t);
pl = conv(xc, 4);
figure
plot(t,xc,'r');
hold on
stem(t,xc);

for t = tmin:dt:tmax;
  xc = cos(2*pi*f*t);
end

xc=zeros(1,20);
n = 0;
for t = tmin:dt:tmax;
  n = n + 1;
  xc(n) = cos(2*pi*f*t);
end

f=90000;
fs = 4*f;
N = 20;
T = 1/f;
w=2*pi*f;
tmax = N*T;
t=linspace(0,tmax,fs*tmax + 1);
x = cos(w*t);
fprintf(' %d\n', x);
figure
plot(x,'r');
hold on
stem(x);

figure
plot(abs(fft(x)))% absolute spectrum of x

% матрица коэффициентов 51 коэффициент
kf = [0 116 0 -214 0 354 0 -540 0 772 0 -1042 0 1339 0 -1645 0 1939 ...
      0 -2199 0 2403 0 -2533 0 2577 0 -2533 0 2403 0 -2199 0 1939 ...
      0 -1645 0 1939 0 -1042 0 772 0 -540 0 354 0 -214 0 116 0];

t=linspace(tmin, tmax+tmin, fs*tmax+1);
% матрицs сигналов
x1 = A*sin(w1*t);

%% инициализация

clc;
clear all;

f1 = 90000; % центральная частота Гц
f2 = 70000; % центральная частота Гц
fs = 360000; % частота дискретизации Гц
A = 2047; % амплитуда сигнала
ts = 1/fs; % период одного дискрета
tm = 50; % количество проходов участка моделирования
tm_test = tm*4; % общее количество проходов моделирования
td = 5; % коэффициент децимации
w1=2*pi*f1; % круговая частота 1
w2=2*pi*f2; % круговая частота 2
tmin = 0; % начальное время (начальная фаза)
tmax = tm*ts; % конечное время
tmax_test = tm_test*ts; % конечное время

% матрица коэффициентов 51 коэффициент (Equiripple 75-85-95-105 кГц)
kf = [0,-834,0,283,0,-202,0,8,0,307,0,-738,0,1260,0,-1836,0,2417, ...
0,-2947,0,3372,0,-3647,0,3742,0,-3647,0,3372,0,-2947,0, ...
2417,0,-1836,0,1260,0,-738,0,307,0,8,0,-202,0,283,0,-834,0];

%lkf =length(kf); % количество коэффициентов
% матрица времени
t=linspace(tmin, tmax_test+tmin, fs*tmax_test+1);
% матрица сигнала 90 кГц
x1 = A*sin(w1*t);
x1(1:60) = 0;
x1(140:201) = 0;
% матрица сигнала 70 кГц
x2 = A*sin(w2*t);
x2(1:60) = 0;
x2(140:201) = 0;

%% моделирование

zr1=zeros(1,200);
zr2=zeros(1,200);
for tm_mod = 1:tm_test

% поэлементное умножение
if tm_mod < 150
z1 = x1(tm_mod:(tm_mod+50)).*kf;
z2 = x2(tm_mod:(tm_mod+50)).*kf;
end
% матрица выходных данных КИХ для 70 кГц и 90 кГц
zr1(tm_mod) = round(sum(z1));
zr2(tm_mod) = round(sum(z2));

end

% максимальное значение
fprintf('Максимальное значение 90 кГц %d\n', max(zr1));
fprintf('Максимальное значение 70 кГц %d\n', max(zr2));
fprintf('\n');

%% графика

% выборки по сигналу x1
figure(1)
plot(x1, 'r');
hold on
stem(x1);

% графика элементов массива zr1
figure(2)
plot(zr1, 'r'), grid;
hold on
stem(zr1);

% графика элементов массива zr2
figure(3)
plot(zr2, 'r');
hold on
stem(zr2);

% absolute spectrum of x1
figure (4)
plot(abs(fft(x1)))