% =====================================================================
% Podstawy techniki cyfrowej, studia niestacjonarne, zjazd 1 (sygnały dyskretne, przetworniki A/C i C/A):
% zadania w Octave (O.1-O.4 obowiązkowe, O.5 dodatkowe)
% =====================================================================
% Ten plik to część zjazdu 1 do zrobienia w Octave (reszta jest w notatniku Jupyter).
% Jak pracować (myCompiler: Octave w przeglądarce, bez konta):
%  1. Otwórz ten plik w zwykłym edytorze tekstu (Notatnik, TextEdit) i skopiuj całą jego zawartość.
%  2. Wejdź na https://[Log in to view URL], usuń przykładowy kod i wklej skopiowany tekst.
%  3. Uzupełniasz tylko linie z komentarzem "% TODO": w miejsce [] wpisujesz właściwe wyrażenie.
%  4. Kliknij Uruchom. Po prawej zobaczysz sprawdzenia [OK], [BŁĄD] albo [DO ZROBIENIA] i wykresy.
%     Dźwięku tam nie posłuchasz: to normalne, do oceny wystarczą sprawdzenia.
%  5. Ostatnią linię wydruku (--- PODSUMOWANIE ...) wklej do notatnika.
%  6. Skopiuj cały kod z powrotem do tego pliku na komputerze i zapisz go. Ten plik oddajesz w Teams
%     razem z notatnikiem, do końca zjazdu (17:00).
% Octave na komputerze też działa: otwórz plik i naciśnij F5.
% Rysunki są tylko trzy, bo myCompiler przerywa program po 15 sekundach, a każdy rysunek trwa tam 1-2 s
% (pozostałe wykresy są w notatniku).
% Najważniejsza różnica względem Pythona: w Octave numeracja zaczyna się od 1, x(1) to pierwsza próbka,
% czyli x[0] ze wzoru.
% =====================================================================
1;  % ta linia mówi Octave, że to skrypt (a funkcje niżej są jego częścią)
% =====================================================================
% FUNKCJE POMOCNICZE (nie trzeba ich czytać ani zmieniać)
% Są w tym samym pliku, żeby wystarczyło wkleić jeden plik do myCompiler (albo otworzyć w Octave).
% =====================================================================

function ptc_start(nazwa)
  % Zeruje liczniki sprawdzeń i ustawia tryb pracy.
  global PTC_WYNIKI PTC_TEST
  PTC_WYNIKI = struct('ok', 0, 'blad', 0, 'brak', 0, 'nazwa', nazwa);
  PTC_TEST = ~isempty(getenv('PTC_TEST'));
  if ~PTC_TEST
    close all;
  end
  printf('\n=== %s ===\n\n', nazwa);
end

function t = gotowe(varargin)
  % Prawda, gdy wszystkie podane zmienne są już uzupełnione (nie są puste).
  t = true;
  for i = 1:numel(varargin)
    if isempty(varargin{i})
      t = false;
    end
  end
end

function sprawdz(opis, czy_gotowe, warunek)
  % Wypisuje [OK], [BŁĄD] albo [DO ZROBIENIA].
  %   opis        - nazwa sprawdzenia, np. '2.1: splot policzony ręcznie'
  %   czy_gotowe  - wynik funkcji gotowe(...): czy zadanie jest uzupełnione
  %   warunek     - funkcja bez argumentów, np. @() isequal(y, [2 5 3 -1 -1])
  global PTC_WYNIKI
  if ~czy_gotowe
    printf('[DO ZROBIENIA]  %s\n', opis);
    PTC_WYNIKI.brak = PTC_WYNIKI.brak + 1;
    return
  end
  try
    ok = all(logical(warunek()));
  catch err
    printf('[BŁĄD]          %s\n                (Octave zgłasza: %s)\n', opis, err.message);
    PTC_WYNIKI.blad = PTC_WYNIKI.blad + 1;
    return
  end
  if ok
    printf('[OK]            %s\n', opis);
    PTC_WYNIKI.ok = PTC_WYNIKI.ok + 1;
  else
    printf('[BŁĄD]          %s\n                (wynik się nie zgadza: przeczytaj jeszcze raz opis zadania w notatniku i komentarz nad zadaniem)\n', opis);
    PTC_WYNIKI.blad = PTC_WYNIKI.blad + 1;
  end
end

function podsumowanie()
  % Wypisuje liczbę zaliczonych sprawdzeń. Tę linię wklejasz przy oddawaniu pracy.
  global PTC_WYNIKI
  w = PTC_WYNIKI;
  printf('\n--- PODSUMOWANIE (%s): OK %d, BŁĄD %d, DO ZROBIENIA %d ---\n', ...
         w.nazwa, w.ok, w.blad, w.brak);
end

function rysunek(szer, wys)
  % Otwiera nowe okno wykresu (opcjonalnie o podanym rozmiarze w pikselach).
  global PTC_TEST
  if nargin < 1, szer = 800; end
  if nargin < 2, wys = 380; end
  if PTC_TEST
    figure('visible', 'off', 'position', [100 100 szer wys]);
  else
    figure('position', [100 100 szer wys]);
  end
end

function zapisz_rysunek(nazwa)
  % W zwykłej pracy tylko odświeża wykres. Przy budowaniu materiałów zapisuje go do pliku PNG.
  global PTC_TEST
  if PTC_TEST
    katalog = getenv('PTC_RYS');
    if isempty(katalog), katalog = '.'; end
    print(gcf, '-dpng', '-r110', fullfile(katalog, [nazwa '.png']));
    close(gcf);
  else
    drawnow;
  end
end

function posluchaj(x, fs, nazwa, normalizuj)
  % Zapisuje dźwięk do pliku nazwa.wav i (jeśli się da) od razu go odtwarza.
  % W myCompiler dźwięku nie słychać (do oceny wystarczą sprawdzenia); w Octave na komputerze plik .wav można odsłuchać.
  % normalizuj = true (domyślnie): głośność wyrównana do maksimum; false: bez zmian (wartości poza [-1, 1] obcięte).
  global PTC_TEST ODTWARZAJ
  if nargin < 4, normalizuj = true; end
  x = double(x(:));
  if normalizuj && max(abs(x)) > 0
    x = 0.9 * x / max(abs(x));
  end
  x = min(max(x, -1), 1);
  audiowrite([nazwa '.wav'], x, fs);
  printf('dźwięk: zapisano plik %s.wav (%.1f s)\n', nazwa, numel(x) / fs);
  if ~PTC_TEST && ~isempty(ODTWARZAJ) && ODTWARZAJ
    try
      sound(x, fs);
    catch
      printf('  (odtwarzanie niedostępne: pobierz plik %s.wav i odsłuchaj go)\n', nazwa);
    end
  end
end

function [f, XdB] = widmo_db(x, fs, okno)
  % Widmo amplitudowe w dB dla częstotliwości 0 ... fs/2, maksimum = 0 dB.
  x = double(x(:));
  if nargin < 3 || isempty(okno)
    okno = ones(size(x));
  end
  N = numel(x);
  X = abs(fft(x .* okno(:)));
  X = X(1:floor(N/2) + 1);
  f = (0:floor(N/2))' * fs / N;
  XdB = 20 * log10(X / max(X) + 1e-12);
end

function spektrogram(x, fs, tytul, fmax)
  % Spektrogram: czas na osi poziomej, częstotliwość na pionowej, kolor = moc w dB.
  if nargin < 3, tytul = ''; end
  if nargin < 4 || isempty(fmax), fmax = fs / 2; end
  x = double(x(:));
  Nf = 512; krok = 128;
  okno = hanning(Nf);
  ile = floor((numel(x) - Nf) / krok) + 1;
  S = zeros(Nf/2 + 1, ile);
  for i = 1:ile
    fragment = x((i-1)*krok + (1:Nf)) .* okno;
    X = abs(fft(fragment));
    S(:, i) = X(1:Nf/2 + 1);
  end
  SdB = 20 * log10(S / max(S(:)) + 1e-12);
  t = ((0:ile-1) * krok + Nf/2) / fs;
  f = (0:Nf/2) * fs / Nf;
  imagesc(t, f, SdB, [-100 0]);
  axis xy; ylim([0 fmax]);
  colormap(jet); c = colorbar; ylabel(c, 'dB');
  xlabel('t [s]'); ylabel('f [Hz]'); title(tytul);
end

function h = projektuj_fir(liczba_wsp, fg, fs, typ, okno, skaluj)
  % Filtr FIR metodą okien (jak fir1 z pakietu signal, ale bez pakietu).
  %   liczba_wsp - liczba współczynników (najlepiej nieparzysta)
  %   fg         - częstotliwość graniczna w Hz; dla pasmowego [f1 f2]
  %   typ        - 'dolny', 'gorny' albo 'pasmowy'
  %   okno       - 'hamming' (domyślnie), 'hanning', 'blackman' albo 'prostokatne'
  %   skaluj     - true (domyślnie): wzmocnienie 1 w środku pasma przepustowego
  if nargin < 4, typ = 'dolny'; end
  if nargin < 5 || isempty(okno), okno = 'hamming'; end
  if nargin < 6, skaluj = true; end
  M = liczba_wsp - 1;
  n = (0:M)' - M/2;
  lp = @(fc) 2 * fc / fs * sinc(2 * fc / fs * n);       % idealny dolnoprzepustowy
  switch typ
    case 'dolny'
      h = lp(fg(1));
    case 'gorny'
      h = (n == 0) - lp(fg(1));
    case 'pasmowy'
      h = lp(fg(2)) - lp(fg(1));
    otherwise
      error('typ filtru: dolny, gorny albo pasmowy');
  end
  switch okno
    case 'hamming',     w = 0.54 - 0.46 * cos(2*pi*(0:M)'/M);
    case 'hanning',     w = 0.5 - 0.5 * cos(2*pi*(0:M)'/M);
    case 'blackman',    w = 0.42 - 0.5 * cos(2*pi*(0:M)'/M) + 0.08 * cos(4*pi*(0:M)'/M);
    case 'prostokatne', w = ones(M+1, 1);
    otherwise, error('okno: hamming, hanning, blackman albo prostokatne');
  end
  h = h .* w;
  if skaluj
    switch typ
      case 'dolny',   f0 = 0;
      case 'gorny',   f0 = fs / 2;
      case 'pasmowy', f0 = mean(fg);
    end
    h = h / abs(sum(h .* exp(-1j * 2*pi*f0/fs * (0:M)')));
  end
  h = h.';
end

function x = samogloski(fs)
  % Sztuczna mowa: samogłoski a, e, i, o, u z syntezy formantowej
  % (ciąg impulsów jak drgania strun głosowych + trzy rezonatory jak jama ustna).
  if nargin < 1, fs = 16000; end
  F = [730 1150 2600; 530 1750 2500; 300 2250 3000; 520 900 2500; 330 750 2350];
  B = [90 110 160];
  f0 = 120; dl = 0.35;
  x = [];
  for i = 1:5
    n = round(dl * fs);
    t = (0:n-1) / fs;
    f_t = f0 * (1.08 - 0.12 * t / dl) * (1 + 0.03 * (i-1));
    faza = cumsum(f_t) / fs;
    y = double(diff([0 floor(faza)]) > 0);
    for k = 1:3
      r = exp(-pi * B(k) / fs); th = 2 * pi * F(i, k) / fs;
      y = filter(1 - r, [1, -2*r*cos(th), r^2], y);
    end
    y = diff([0 y]);
    y = y .* sin(pi * (0:n-1) / n) .^ 0.5;
    x = [x, y, zeros(1, round(0.08 * fs))];
  end
  x = 0.8 * x / max(abs(x));
end

function x = melodia(fs, harmoniczne)
  % Melodia „Wlazł kotek na płotek”: tony z harmonicznymi o amplitudach 1/k.
  if nargin < 1, fs = 16000; end
  if nargin < 2, harmoniczne = 12; end
  nuty = [392.0 329.6 329.6 349.2 293.7 293.7 261.6 329.6 392.0];
  n = round(0.3 * fs);
  t = (0:n-1) / fs;
  obw = min(t / 0.01, 1) .* exp(-4 * t);
  x = [];
  for f = nuty
    ton = zeros(1, n);
    for k = 1:harmoniczne
      if k * f < fs / 2
        ton = ton + sin(2*pi*k*f*t) / k;
      end
    end
    x = [x, ton .* obw];
  end
  x = 0.8 * x / max(abs(x));
end

function x = nuta(fs, f0)
  % Jeden dźwięk instrumentu: podstawowa najsilniejsza, wyższe harmoniczne słabsze i szybciej gasnące.
  if nargin < 1, fs = 16000; end
  if nargin < 2, f0 = 220; end
  t = (0:round(1.5 * fs)-1) / fs;
  a = [1 0.5 0.35 0.25 0.15 0.1];
  x = zeros(size(t));
  for k = 1:6
    x = x + a(k) * exp(-k * 1.5 * t) .* sin(2*pi*k*f0*t);
  end
  x = x .* min(t / 0.005, 1);
  x = 0.8 * x / max(abs(x));
end

% =====================================================================
% KONIEC FUNKCJI POMOCNICZYCH. Właściwe zadania są niżej.
% =====================================================================

global ODTWARZAJ
ODTWARZAJ = 0;   % 0 = dźwięki tylko zapisywane do plików .wav (działa wszędzie);
                 % 1 = odtwarzaj od razu (tylko Octave na komputerze; w przeglądarce daje błędy ALSA)
ptc_start('Zjazd 1: zadania w Octave');

%% ---------------------------------------------------------------------
%% Zadanie O.1: delta, skok, impuls prostokątny i ton 440 Hz
%% ---------------------------------------------------------------------
% To samo co zadania 1.1, 1.2 i 3.1 w notatniku, ale w Octave. Pięć linii z TODO O.1.
% Funkcję w jednej linii zapisuje się jako funkcję anonimową:  nazwa = @(n) wyrażenie;
% na przykład kwadrat = @(n) n.^2;
% Porównanie n == 5 daje 1 tam, gdzie warunek jest spełniony, i 0 gdzie indziej; double(...) zamienia
% wynik porównania na liczby. Przykład: piatka = @(n) double(n == 5);
%  - delta: jedynka tylko dla n = 0;  skok: jedynki dla n >= 0 (porównanie "większe lub równe" to >=);
%  - p: impuls prostokątny u[n] - u[n-5]: od skoku odejmij skok przesunięty o 5, czyli skok(n - 5)
%    (tu bez @(n): p to już liczby, a nie funkcja);
%  - t: chwile próbkowania 0, 1/fs, 2/fs, ...: zapis (0:3) / 4 daje 0 0.25 0.5 0.75 (nawias jest ważny);
%  - x: ton A cos(2 pi f t); cos i pi bez np., cos działa na całym wektorze naraz.
n = -10:10;
delta = @(n) double(n == 0);  % TODO O.1: delta Kroneckera: 1 dla n = 0, w pozostałych miejscach 0
skok = @(n) double(n >= 0);   % TODO O.1: skok jednostkowy: 1 dla n >= 0, 0 dla n < 0
p = [];
if gotowe(delta, skok)
p = skok(n) - skok(n - 5);  % TODO O.1: impuls prostokątny u[n] - u[n-5]
end
if gotowe(delta, skok, p)
  rysunek(900, 260);
  subplot(1, 3, 1); stem(n, delta(n), 'filled'); title('\delta[n]'); xlabel('n'); ylim([-0.2 1.2]);
  subplot(1, 3, 2); stem(n, skok(n), 'filled'); title('u[n]'); xlabel('n'); ylim([-0.2 1.2]);
  subplot(1, 3, 3); stem(n, p, 'filled'); title('u[n] - u[n-5]'); xlabel('n'); ylim([-0.2 1.2]);
  zapisz_rysunek('Z1_O1_sygnaly');
end
sprawdz('O.1: delta', gotowe(delta), @() isequal(delta(-2:2), [0 0 1 0 0]));
sprawdz('O.1: skok', gotowe(skok), @() isequal(skok(-2:2), [0 0 1 1 1]));
sprawdz('O.1: impuls ma 5 jedynek, od n = 0 do n = 4', gotowe(p), ...
        @() sum(p) == 5 && p(n == 0) == 1 && p(n == 4) == 1 && p(n == 5) == 0);

fs = 8000; f = 440; A = 0.5;                  % ton 440 Hz, amplituda 0,5, 8000 próbek na sekundę, 1 s
t = (0:7999) / fs;            % TODO O.1: chwile próbkowania 0, 1/fs, 2/fs, ... (8000 liczb)
x = [];
if gotowe(t)
x = A * cos(2 * pi * f * t); % TODO O.1: sinusoida 440 Hz
end
if gotowe(t, x)
  posluchaj(x, fs, 'Z1_ton_440Hz');
end
sprawdz('O.1: 8000 próbek co 1/8000 s', gotowe(t), @() numel(t) == 8000 && abs(t(2) - 1/8000) < 1e-12);
sprawdz('O.1: wartości sinusoidy', gotowe(x), ...
        @() abs(x(1) - 0.5) < 1e-12 && abs(x(11) - 0.5 * cos(2*pi*440*10/8000)) < 1e-12);

%% ---------------------------------------------------------------------
%% Zadanie O.2: te same próbki dwóch częstotliwości i zespolona sinusoida (funkcja kołowa)
%% ---------------------------------------------------------------------
% To samo co zadania 3.2, 4.1 i 4.2 w notatniku. Cztery linie z TODO O.2.
%  - x1, x2: cosinusy o częstotliwości 0,3 i 1,3 cyklu na próbkę dla nn = 0, 1, ..., 19. Wzór jak w O.1,
%    tylko zamiast f * t jest od razu częstotliwość razy nn. Przykład dla 0,1 cyklu na próbkę:
%    cos(2 * pi * 0.1 * nn). Część dziesiętną oddziela kropka.
%  - z: zespolona sinusoida e^(j w0 k): funkcja exp(...), a w nawiasie jednostka urojona 1j (cyfra jeden
%    i litera j) razy w0 razy kk. Przykład na innej liczbie: exp(1j * pi / 2) daje (prawie dokładnie) 1j.
%  - c: cosinus z dwóch wskazów (wzór Eulera): suma z i liczby sprzężonej conj(z), podzielona przez 2.
%    Przykład na liczbach: (3 + 4j) i jej sprzężenie (3 - 4j) dają razem 6, a połowa to 3.
%    W Octave numeracja jest od 1: z(1) to z[0] ze wzoru, a z(5) to z[4].
nn = 0:19;
x1 = cos(2 * pi * 0.3 * nn);  % TODO O.2: cosinus 0,3 cyklu na próbkę
x2 = cos(2 * pi * 1.3 * nn);  % TODO O.2: cosinus 1,3 cyklu na próbkę
kk = 0:31;
w0 = 2 * pi / 16;
z = exp(1j * w0 * kk);  % TODO O.2: zespolona sinusoida, wskaz obracający się o w0
c = [];
if gotowe(z)
c = (z + conj(z)) / 2;  % TODO O.2: cosinus z sumy dwóch wskazów
end
if gotowe(z)
  rysunek(500, 420);
  fi = linspace(0, 2 * pi, 200);
  plot(cos(fi), sin(fi), 'color', [0.7 0.7 0.7]); hold on;
  plot(real(z), imag(z), 'bo', 'markerfacecolor', 'b'); hold off; axis equal;
  xlabel('część rzeczywista'); ylabel('część urojona'); title('z[k] = e^{j w_0 k}: punkty na okręgu');
  zapisz_rysunek('Z1_O2_wskaz');
end
sprawdz('O.2: obie sinusoidy dają te same 20 próbek', gotowe(x1, x2), ...
        @() numel(x1) == 20 && max(abs(x1 - x2)) < 1e-9);
sprawdz('O.2: |z| = 1, po 4 próbkach z = j, po 16 próbkach z = 1', gotowe(z), ...
        @() numel(z) == 32 && max(abs(abs(z) - 1)) < 1e-12 && abs(z(5) - 1j) < 1e-12 && abs(z(17) - 1) < 1e-12);
sprawdz('O.2: (z + conj(z)) / 2 = cos(w0 k)', gotowe(c), ...
        @() max(abs(c - cos(w0 * kk))) < 1e-12);

%% ---------------------------------------------------------------------
%% Zadanie O.3: przetworniki C/A i A/C jako funkcje, Arduino mierzy 3,3 V
%% ---------------------------------------------------------------------
% To samo co zadania 6.2 i 6.3 w notatniku. Pięć linii z TODO O.3: po jednej w funkcjach ca i ac
% i trzy niżej.
%  - krok przetwornika to Uref podzielone przez 2 do potęgi B (potęga w Octave to ^): dla 4 V i 2 bitów
%    4 / 2^2 daje 1 V;
%  - ca (kod -> napięcie): kod razy krok. Przykład: kody [0 1 3] przy kroku 1 V dają 0, 1 i 3 V;
%  - ac (napięcie -> kod): napięcie podzielone przez krok i zaokrąglone W DÓŁ funkcją floor(...)
%    (krok weź w nawias, bo sam jest dzieleniem). Przykład: floor(2.7 / 1) daje 2. Obcięcie do zakresu
%    0 ... 2^B - 1 jest gotowe w następnej linii;
%  - kod_33: kod Arduino (10 bitów, 5 V) dla 3,3 V, czyli ac(...) z tymi trzema liczbami;
%    U_33: napięcie odtworzone z tego kodu funkcją ca(...); blad_mV: prawdziwe napięcie minus odtworzone,
%    razy 1000. Przykład zapisu dla innego napięcia: ac(1.2, 10, 5.0).
function U = ca(kod, B, Uref)
  U = kod * (Uref / 2^B);  % TODO O.3: kod -> napięcie (kod razy krok)
end

function k = ac(U, B, Uref)
  k = floor(U / (Uref / 2^B));  % TODO O.3: napięcie -> numer poziomu
  k = min(max(k, 0), 2^B - 1);  
end
kod_33 = ac(3.3, 10, 5.0);  % TODO O.3: kod Arduino dla 3,3 V (10 bitów, 5 V)
U_33 = [];
blad_mV = [];
if gotowe(kod_33)
 U_33 = ca(kod_33, 10, 5.0); % TODO O.3: napięcie odtworzone z kodu
end
if gotowe(U_33)
  blad_mV = (3.3 - U_33) * 1000; % TODO O.3: błąd w miliwoltach
end
if gotowe(ca(5, 3, 8.0), ac(2.6, 3, 8.0))
  printf('Arduino: 2,5 V -> %d   kod 1023 -> %.4f V\n', ac(2.5, 10, 5.0), ca(1023, 10, 5.0));
end
if gotowe(blad_mV)
  printf('3,3 V: kod %d, napięcie odtworzone %.4f V, błąd %.2f mV, krok %.2f mV\n', kod_33, U_33, blad_mV, 5 / 1024 * 1000);
end
sprawdz('O.3: C/A, kod 5 przy 3 bitach i 8 V daje 5 V', gotowe(ca(5, 3, 8.0)), @() abs(ca(5, 3, 8.0) - 5) < 1e-12);
sprawdz('O.3: A/C, 2,6 V -> 2; 9 V -> 7 (nasycenie); -1 V -> 0', gotowe(ac(2.6, 3, 8.0)), ...
        @() isequal(ac([2.6 9.0 -1.0], 3, 8.0), [2 7 0]));
sprawdz('O.3: Arduino, 2,5 V -> 512', gotowe(ac(2.5, 10, 5.0)), @() ac(2.5, 10, 5.0) == 512);
sprawdz('O.3: 3,3 V -> kod 675, napięcie odtworzone i błąd', gotowe(kod_33, U_33, blad_mV), ...
        @() kod_33 == 675 && abs(U_33 - 675 * 5 / 1024) < 1e-12 && abs(blad_mV - (3.3 - 675 * 5 / 1024) * 1000) < 1e-9);

%% ---------------------------------------------------------------------
%% Zadanie O.4: kwantyzator B-bitowy i stosunek sygnału do szumu kwantyzacji (SQNR)
%% ---------------------------------------------------------------------
% To samo co zadania 7.1 i 7.2 w notatniku. Trzy linie z TODO O.4: dwie w funkcji kwantyzuj, jedna sqnr.
% Funkcja zapisana w skrypcie wygląda tak: function wynik = nazwa(argumenty) ... end.
%  - d: krok to długość zakresu dźwięku (od -1 do 1, czyli 2) podzielona przez liczbę poziomów,
%    czyli przez 2 do potęgi B. Potęga w Octave to ^ (nie ** jak w Pythonie): 2^3 daje 8,
%    więc przy 3 bitach d wychodzi 0,25.
%  - xq: trzy ruchy jak w zadaniu 7.1: x dzielisz przez krok d, zaokrąglasz funkcją round(...)
%    (bez np.) i mnożysz z powrotem przez d. Przykład: krok 0,25 i wartość 0,6 dają 0,5.
%    Wartości spoza zakresu obcina do [-1, 1 - d] gotowa ostatnia linia funkcji.
%  - sqnr: funkcja w jednej linii, która zaczyna się od @(x, xq), a po niej jest wyrażenie
%    (przykład innej takiej funkcji: podwoj = @(x) 2 * x;). Na górze ułamka suma kwadratów dźwięku,
%    na dole suma kwadratów błędu (różnicy x - xq, w nawiasie). Górę dzielisz przez dół, z wyniku
%    logarytm dziesiętny log10(...), całość razy 10. Kwadrat każdego elementu to .^ 2 (z kropką),
%    suma to sum(...). Przykład: dla sum 9 i 0,25 wychodzi około 15,6 dB.
function xq = kwantyzuj(x, B)
  d = 2 / 2^B;  % TODO O.4: krok kwantyzacji
  xq = round(x / d) * d;  % TODO O.4: zaokrąglenie x do najbliższego poziomu
  xq = min(max(xq, -1), 1 - d);
end

fs = 8000;
t = (0:31) / fs;                              % 32 próbki = 4 ms
xk = 0.95 * sin(2 * pi * 500 * t);
sprawdz('O.4: 0,3 przy 2 bitach -> 0,5', gotowe(kwantyzuj(xk, 2)), @() abs(kwantyzuj(0.3, 2) - 0.5) < 1e-12);
sprawdz('O.4: 1,0 przy 3 bitach -> 0,75 (nasycenie)', gotowe(kwantyzuj(xk, 2)), @() abs(kwantyzuj(1.0, 3) - 0.75) < 1e-12);
sprawdz('O.4: poza nasyceniem błąd nie większy niż pół kroku', gotowe(kwantyzuj(xk, 2)), ...
        @() max(abs(xk(abs(xk) < 0.8) - kwantyzuj(xk(abs(xk) < 0.8), 4))) <= 2 / 2^4 / 2 + 1e-12);

sqnr = @(x, xq) 10 * log10(sum(x .^ 2) / sum((x - xq) .^ 2));  % TODO O.4: SQNR w dB

fs = 48000;
t = (0:fs-1) / fs;                            % 1 s
sinus = sin(2 * pi * 997 * t);                % 997 Hz: błąd kwantyzacji nie powtarza się okresowo
bity = 2:16;
teoria = 6.02 * bity + 1.76;
zmierzone = [];
if gotowe(kwantyzuj(sinus, 2), sqnr)
  zmierzone = zeros(size(bity));
  for i = 1:numel(bity)
    B = bity(i);
    xa = (1 - 2 / 2^B) * sinus;               % amplituda 1 - d: największa bez nasycenia
    zmierzone(i) = sqnr(xa, kwantyzuj(xa, B));
    printf('B = %2d:  zmierzone %6.2f dB,  teoria %6.2f dB\n', B, zmierzone(i), teoria(i));
  end
  rysunek(700, 380);
  plot(bity, teoria, 'b-', 'linewidth', 1.5); hold on;
  plot(bity, zmierzone, 'ro', 'markerfacecolor', 'r'); hold off; grid on;
  xlabel('B [bity]'); ylabel('SQNR [dB]'); legend('teoria 6,02 B + 1,76', 'zmierzone', 'location', 'northwest');
  title('każdy bit to około 6 dB');
  zapisz_rysunek('Z1_O4_sqnr');
end
sprawdz('O.4: SQNR dla x = [1 1], xq = [1 0] to 3,01 dB', gotowe(sqnr), @() abs(sqnr([1 1], [1 0]) - 10 * log10(2)) < 1e-12);
sprawdz('O.4: SQNR zgodny z teorią (do 0,2 dB) od 6 bitów w górę', gotowe(zmierzone), ...
        @() max(abs(zmierzone(bity >= 6) - teoria(bity >= 6))) < 0.2);

%% ---------------------------------------------------------------------
%% Zadanie O.5 (dodatkowe, +0,5 pkt): liczby int16: nasycenie w Octave, zawijanie w prostym procesorze
%% ---------------------------------------------------------------------
% To zadanie jest DODATKOWE: daje +0,5 pkt (nadrabia braki, razem najwyżej 20 pkt). Jeśli go nie robisz,
% jego 3 sprawdzenia zostaną [DO ZROBIENIA], a w PODSUMOWANIU będzie OK 17, DO ZROBIENIA 3: to jest w porządku.
% W pliku WAV i w rejestrach procesora próbki to liczby całkowite int16: od -32768 do 32767.
% W notatniku (część 9) NumPy "zawijał" wynik spoza zakresu: 20000 * 2 dawało -25536.
% Krok 1 (a): sprawdź, co robi Octave. Policz 20000 razy 2, ale liczbę 20000 zapisz jako typ int16
%         (funkcja int16(...)). Wynik zobaczysz w wydruku.
% Krok 2 (zawin): zasymuluj zawijanie prostego procesora jako funkcję w jednej linii, która zaczyna się
%         od @(y). Trzy ruchy: dodajesz do y liczbę 32768, z wyniku bierzesz resztę z dzielenia przez 65536
%         funkcją mod(...) (w nawiasie po kolei: co dzielisz i przez co; mod(7, 5) daje 2), na końcu
%         odejmujesz 32768. Dla y = 40000 ma wyjść -25536, dla y = 1000 zostaje 1000.
% Krok 3 (y_zaw): wzmocnij melodię 4 razy "jak prosty procesor": m16 to melodia jako liczby int16.
%         Najpierw zamień m16 na zwykłe liczby funkcją double(...), bo bez tego Octave nasyci wynik,
%         zanim zadziała zawijanie. Potem wzmacniasz ją 4 razy, przepuszczasz przez swoją funkcję
%         zawin(...) i dzielisz przez 32768, żeby wrócić do zakresu od -1 do 1.
a = int16(20000) * 2;  % TODO O.5: 20000 jako int16 razy 2
zawin = @(y) mod(y + 32768, 65536) - 32768;  % TODO O.5: zawijanie do zakresu -32768 ... 32767
fs = 16000;
m = melodia(fs);
m16 = int16(round(m * 32767));                % melodia jako liczby int16 (jak w pliku WAV)
y_oct = double(int16(4) * m16) / 32767;       % wzmocnienie x4 policzone przez Octave na int16 (gotowe)
y_zaw = [];
if gotowe(zawin)
y_zaw = zawin(4 * double(m16)) / 32768;  % TODO O.5: wzmocnienie x4 z zawijaniem
end
if gotowe(a)
  printf('Octave: int16(20000) * 2 = %d\n', a);
end
if gotowe(y_zaw)
  printf('melodia x4 w Octave na int16: najmniejsza %.4f, największa %.4f (nasycenie)\n', min(y_oct), max(y_oct));
  printf('melodia x4 z zawijaniem: najmniejsza %.4f, największa %.4f\n', min(y_zaw), max(y_zaw));
  posluchaj(y_oct, fs, 'Z1_nasycenie', false);
  posluchaj(y_zaw, fs, 'Z1_przepelnienie', false);
end
sprawdz('O.5: Octave nasyca: int16(20000) * 2 = 32767', gotowe(a), @() isa(a, 'int16') && a == 32767);
sprawdz('O.5: zawijanie: 40000 -> -25536, -40000 -> 25536, 1000 -> 1000', gotowe(zawin), ...
        @() isequal(zawin([40000 -40000 1000]), [-25536 25536 1000]));
sprawdz('O.5: melodia x4 z zawijaniem (skoki o całą skalę, małe wartości bez zmian)', gotowe(y_zaw), ...
        @() numel(y_zaw) == numel(m) && max(abs(y_zaw)) <= 1 && any(abs(diff(y_zaw)) > 1) && ...
            max(abs(y_zaw(abs(m) < 0.2) - 4 * double(m16(abs(m) < 0.2)) / 32768)) < 1e-12);

podsumowanie();