clc
 
close 
clear
% узкополосная фильтрация
x0 = 1.5;
w0 = 2*pi*60;
t = transpose(0:0.001:0.8);
N = 4096;
 
n = length(t);
Y = zeros(n,1);
for ii = 1:n
    if t(ii) >= 0 && t(ii) <= 0.4
        Y(ii) = 0;
    elseif t(ii) > 0.4  && t(ii) <= 0.8 
        Y(ii) = x0*cos(w0*t(ii));
    end
end
figure (1)
plot (t,Y)
 
% проводим дискретизацию по времени
n1 = transpose(0:N-1);
tt = n1.*(t(end)/N);
Y1 = zeros(n,1);
for ii = 1:length(n1)
    if tt(ii) >= 0 && tt(ii) <= 0.4
        Y1(ii) = 0;
    elseif tt(ii) > 0.4  && tt(ii) <= 0.8 
        Y1(ii) = x0*cos(w0*tt(ii));
    end
end
figure (2)
plot(tt,Y1)
 
 
% накладываем гауссовский белый шум
sigma = 2;
sz = size(n1);
Noize = normrnd(0,sigma,sz(1),sz(2));
Z = Y1+Noize;
Z1 = Z-Y1;
figure (3)
plot(tt,Z,'k-')
 
 
% подвергаем преобразованию Фурье
%L = [1:N-1]';
%nfft = length(L);
ns = (length(tt));
ks=1/tt(2);
YY = (abs(fft(Y1))/(ns));
YY=YY(1:ns/2+1);
YY(2:end-1) = 2*YY(2:end-1);
ZZ = (abs(fft(Z))/(ns));
ZZ=ZZ(1:ns/2+1);
ZZ(2:end-1) = 2*ZZ(2:end-1);
k=ks*(0:ns/2)/ns;
figure (4)
grid on;
plot(k,YY/2);
hold on
stem(k,ZZ/2);
hold on;
xlim([0 120])
 
% обнуляем те элменты к спектральной функции...
%  к, где спектр полезного сигнала  не  сильно отличается от нуля, и
%  совершаем обратное преобразование Фурье
YY1 = zeros(1,length(ns));
for i = 1:length(ns)
    if ns(i) >= 40 && ns(i) <= 80
        YY1(i) = YY;
    elseif ns(i) < 40  && ns(i) > 80 
        YY1(i) = 0;
    end
end
% обатное преобразование
y = ifft(YY1,ns);
figure (5)
plot(tt,Z,'k-',tt,y+6,'r--')  %,tt,y+6,'r--')
 
 
%строим две АКФ
ZZ1 = Z-Y1; % чистый шум
 
kk = 1:ns/2;
RR = zeros(1,(ns/2));
for k = 1:ns/2
    for i = 1:(N-ns/2-i)
    RR(k) = Z(i)*Z(i+k);
    end
end
% RR = (1/(N-k))*RR;
 
 
RR1 = zeros(1,(ns/2));
for k = 1:ns/2
    for i = 1:(N-ns/2-i)
    RR1(k) = ZZ1(i)*ZZ1(i+k);
    end
end
% RR1 = (1/(N-k))*RR1;
plot (k,RR1,k,RR1);
 
 
 
% Corr_auto1 = autocorr(Z);
% Corr_auto2 = autocorr(Z1);
% figure (6)
% r1 = autocorr(Z,'NumLags',300,'NumSTD',2);
% r2 = autocorr(Z1,'NumLags',300,'NumSTD',2);
% plot (1:length(r1),r1,1:length(r2),r2)
% ylim([-1 1]);
 
% plot (length(Corr_auto1),Corr_auto1,length(Corr_auto2),Corr_auto2);
 
 
 
% clear
% clc
TMAX=1;
N=4096;
dt=TMAX/N;
tn=0:dt:TMAX;
A=1;
% Сигнал
% sig=A*(tn>= 0.5)-(A*(tn >= 1));  %+(-A*(tn >= 1));
sig=A*(tn>=0.5).*(exp(-(5*(tn-0.5).*(tn>=0.5))));
% Сигнал + помеха
x=sig+1*randn(N+1,1)'; % сигнал с шумом
x1 = x-sig; % шум без полезного сигнала
L=length(sig); %длина паттерна
N=length(x); %длина входного сигнала
y=zeros(1,N);
y1=zeros(1,N);
for i=1:N
    for p=1:L
        if (i-p)>0
        %y(i)=y(i)+x(i-p)*sig(N-p)*dt;
%         y(i)=y(i-p)+x(p)*dt;
%         y1(i)=y(i-p)+x1(p)*dt;
        y(p)=y(i-p)+x(p)*dt;
        y1(p)=y(i-p)+x1(p)*dt;
        end
    end
end
figure(7)
% вывод результата
subplot (4,1,1); plot(tn,sig);grid;title('signal') % шаблон сигнала
subplot (4,1,2); plot(tn,x);grid;title('signal+noise') % реальный сигнал
subplot (4,1,3); plot(tn,y);grid;title('filter output') % выход фильтр
subplot (4,1,4); plot(tn,y1);grid;title('filter output without noise') % выход фильтра без учета шума