clear;
N = 4000; %ilosc problek
Nf = 4096;  %Powinnno byc wieksze od ilosci probek
fs = 1000; % frequency sample
t = 0: 1/fs : (N-1)/fs; % okres
A1 = 6;
A2 = 2;% Amplituda
f1 =100;% 100 Hz = częstotliwość
f2 = 400;
e = randn([1,N]); %szum gaussa

x1 = A1 * sin(2*pi*f1*t);
x2 = A2 * sin(2*pi*f2*t);
x3 = x1+x2+e; % Dodanie sygnałów do siebie z zakłóceniem
%fft - fast fourier transfrotm

w = abs(fft(x3,Nf)); %w = amplituda NF musi byc 2^x; Widmo amplitudowe !!WAZNE
N21 = Nf/2 +1; %połowa próblek. +1 bo 0 określa stałą widma  !!!WAZNE
%f= linspace(0,fs,Nf); %Generuje wektor próbek
f = linspace(0,fs/2,N21); %wektor połowy próbek !!WAZNE
figure(1);
plot(t,x3);
figure(2);
plot(f,w(1:N21)); %Wybieramy próbki od 1 do wartości N21
xlabel('Czas');
ylabel('x(t)');
%chirp - Sygnal sinusoidalny o modulowanej czestotliwosci
%arg wektor czasu, czestotliwosc, koniec okna czasowego, f2 dla konca okna czasowego
figure(3);
%sublot, a - obrazkow w wierszach, b - liczba kolumn, c - index obrazka do
%ktorego sie odowlujemy
subplot(3,2,1);
plot(t,x3);
ch = chirp(t,f1,(N-1)/fs,f2);
w1 = abs(fft(ch,Nf));
subplot(3,2,2)
plot(f,w1(1:N21));

% FILTRACJA - DOlna przepustowość
% Założenia: probka {xk, xk-1, xk-2,...}
%Wzór [a0*xk,a1*xk-1,...] ten wektor = fir1;
M=201; % rząd - dlugosc filtru
%drug arg - czestotliwosc graniczna
h =  fir1(M-1,250/(fs/2));
th = 0 : 1/fs:(M-1)/fs;
subplot(3,2,3);
plot(th,h); %h(f)

w2 = abs(fft(h,Nf));
subplot(3,2,4);
plot(f,w2(1:N21)); %widmo h(f)

s = filter(h,1,ch);
subplot(3,2,5);
plot(t,s); %s(t)

w3 = abs(fft(s,Nf));
subplot(3,2,6);
plot(f,w3(1:N21)); %widmo s(t)