function [ ] = zad4()
fs = 8000;

b = oblicz_b(fs, 1000);
[h1, t1] = freqz(b, oblicz_a(0, fs, 1000));
[h2, t2] = freqz(b, oblicz_a(0.5, fs, 1000));
[h3, t3] = freqz(b, oblicz_a(0.9, fs, 1000));
[h4, t4] = freqz(b, oblicz_a(1, fs, 1000));
[h5, t5] = freqz(b, oblicz_a(1.1, fs, 1000));
[h6, t6] = freqz(b, oblicz_a(2, fs, 1000));
[h7, t7] = freqz(b, oblicz_a(3, fs, 1000));

figure();
plot(oblicz_t(t1, fs), oblicz_mag(h1), oblicz_t(t2, fs), oblicz_mag(h2), oblicz_t(t3, fs), oblicz_mag(h3), oblicz_t(t4, fs), oblicz_mag(h4), oblicz_t(t5, fs), oblicz_mag(h5), oblicz_t(t6, fs), oblicz_mag(h6), oblicz_t(t7, fs), oblicz_mag(h7));
title('Charakterystyka amplitudowa, f0 = 1000 Hz');
ylabel('Amplituda [dB]');
xlabel('Częstotliwość [Hz]');
ylim([-75 15]);
legend('r = 0', 'r = 0.5', 'r = 0.9', 'r = 1', 'r = 1.1', 'r = 2', 'r = 3');
grid on;


b = oblicz_b(fs, 2000);
[h1, t1] = freqz(b, oblicz_a(0, fs, 2000));
[h2, t2] = freqz(b, oblicz_a(0.5, fs, 2000));
[h3, t3] = freqz(b, oblicz_a(0.9, fs, 2000));
[h4, t4] = freqz(b, oblicz_a(1, fs, 2000));
[h5, t5] = freqz(b, oblicz_a(1.1, fs, 2000));
[h6, t6] = freqz(b, oblicz_a(2, fs, 2000));
[h7, t7] = freqz(b, oblicz_a(3, fs, 2000));

figure();
plot(oblicz_t(t1, fs), oblicz_mag(h1), oblicz_t(t2, fs), oblicz_mag(h2), oblicz_t(t3, fs), oblicz_mag(h3), oblicz_t(t4, fs), oblicz_mag(h4), oblicz_t(t5, fs), oblicz_mag(h5), oblicz_t(t6, fs), oblicz_mag(h6), oblicz_t(t7, fs), oblicz_mag(h7));
title('Charakterystyka amplitudowa, f0 = 2000 Hz');
ylabel('Amplituda [dB]');
xlabel('Częstotliwość [Hz]');
ylim([-75 15]);
legend('r = 0', 'r = 0.5', 'r = 0.9', 'r = 1', 'r = 1.1', 'r = 2', 'r = 3');
grid on;


b = oblicz_b(fs, 3000);
[h1, t1] = freqz(b, oblicz_a(0, fs, 3000));
[h2, t2] = freqz(b, oblicz_a(0.5, fs, 3000));
[h3, t3] = freqz(b, oblicz_a(0.9, fs, 3000));
[h4, t4] = freqz(b, oblicz_a(1, fs, 3000));
[h5, t5] = freqz(b, oblicz_a(1.1, fs, 3000));
[h6, t6] = freqz(b, oblicz_a(2, fs, 3000));
[h7, t7] = freqz(b, oblicz_a(3, fs, 3000));

figure();
plot(oblicz_t(t1, fs), oblicz_mag(h1), oblicz_t(t2, fs), oblicz_mag(h2), oblicz_t(t3, fs), oblicz_mag(h3), oblicz_t(t4, fs), oblicz_mag(h4), oblicz_t(t5, fs), oblicz_mag(h5), oblicz_t(t6, fs), oblicz_mag(h6), oblicz_t(t7, fs), oblicz_mag(h7));
title('Charakterystyka amplitudowa, f0 = 3000 Hz');
xlabel('Częstotliwość [Hz]');
ylabel('Amplituda [dB]');
ylim([-75 15]);
legend('r = 0', 'r = 0.5', 'r = 0.9', 'r = 1', 'r = 1.1', 'r = 2', 'r = 3');
grid on;

gauss = normrnd(0, sqrt(0.5), 1, fs);

y = filter(oblicz_b(fs, 1000), oblicz_a(0.9, fs, 1000), gauss);
figure();
spectrogram(y, window(@hann, 256), 0, 256, fs, 'yaxis');
title(['Spektrogram filtrowanego szumu, f0 = 1000' ', r = 0.9']);
ylabel(colorbar, 'Moc/częstotliwość [dB/Hz]');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');

y = filter(oblicz_b(fs, 2000), oblicz_a(0.9, fs, 2000), gauss);
figure();
spectrogram(y, window(@hann, 256), 0, 256, fs, 'yaxis');
title(['Spektrogram filtrowanego szumu, f0 = 2000' ', r = 0.9']);
ylabel(colorbar, 'Moc/częstotliwość [dB/Hz]');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');


y = filter(oblicz_b(fs, 3000), oblicz_a(0.9, fs, 3000), gauss);
figure();
spectrogram(y, window(@hann, 256), 0, 256, fs, 'yaxis');
title(['Spektrogram filtrowanego szumu, f0 = 3000' ', r = 0.9']);
ylabel(colorbar, 'Moc/częstotliwość [dB/Hz]');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');

y = filter(oblicz_b(fs, 2000), oblicz_a(0, fs, 2000), gauss);
figure();
spectrogram(y, window(@hann, 256), 0, 256, fs, 'yaxis');
title(['Spektrogram filtrowanego szumu, f0 = 2000' ', r = 0']);
ylabel(colorbar, 'Moc/częstotliwość [dB/Hz]');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');

y = filter(oblicz_b(fs, 2000), oblicz_a(0.5, fs, 2000), gauss);
figure();
spectrogram(y, window(@hann, 256), 0, 256, fs, 'yaxis');
title(['Spektrogram filtrowanego szumu, f0 = 2000' ', r = 0.5']);
ylabel(colorbar, 'Moc/częstotliwość [dB/Hz]');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');

y = filter(oblicz_b(fs, 2000), oblicz_a(1, fs, 2000), gauss);
figure();
spectrogram(y, window(@hann, 256), 0, 256, fs, 'yaxis');
title(['Spektrogram filtrowanego szumu, f0 = 2000' ', r = 1']);
ylabel(colorbar, 'Moc/częstotliwość [dB/Hz]');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');

y = filter(oblicz_b(fs, 2000), oblicz_a(1.1, fs, 2000), gauss);
figure();
spectrogram(y, window(@hann, 256), 0, 256, fs, 'yaxis');
title(['Spektrogram filtrowanego szumu, f0 = 2000' ', r = 1.1']);
ylabel(colorbar, 'Moc/częstotliwość [dB/Hz]');
xlabel('Czas [s]');
ylabel('Częstotliwość [Hz]');
end

function [ t ] = oblicz_t(t, fs)
t = t*fs/(2*pi);
end
function [ h ] = oblicz_mag(h_in)
h = mag2db(abs(h_in));
end
function [ b ] = oblicz_b(fs, f0)
b = [1, -2*cos(2*pi*f0/fs), 1];
end
function [ a ] = oblicz_a(r, fs, f0)
a = [1, -2*r*cos(2*pi*f0/fs), r*r];
end