2 d complete

Old/1_fft_basics/1.py

@@ -0,0 +1,49 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.fft import fft, ifft, fftfreq
+p = True
+N = 100 
+x, dx = np.linspace(0, 2*np.pi, N, endpoint=False, retstep=True)
+f_x = 1 + 2*np.sin(x) + 10*np.sin(5*x) + 3*np.cos(5*x)
+F_k = fft(f_x)
+if p:
+    print("dx:", dx)
+    print("x:", x)
+    print("f(x):", f_x)
+    print("F_k:", F_k)
+
+f_x_reconstructed = ifft(F_k)
+
+# --- VYKRESLENÍ VÝSLEDKŮ PRO KONTROLU ---
+plt.figure(figsize=(12, 8))
+
+# Graf 1: Původní a složená funkce (měly by se překrývat)
+plt.subplot(2, 1, 1)
+plt.plot(x, f_x, label="Původní f(x)", linewidth=3)
+# ifft vrací komplexní čísla, pro graf vezmeme reálnou část (.real)
+plt.plot(x, f_x_reconstructed.real, '--', label="Složená (iFFT)", color='red')
+plt.title("Původní signál vs. Rekonstruovaný signál (Úspěšný test!)")
+plt.legend()
+plt.grid()
+
+# Graf 2: Zobrazení frekvenčního spektra (co vlastně FFT našla)
+plt.subplot(2, 1, 2)
+# Spočítáme frekvence pro osu X (zajímají nás jen kladné frekvence, proto úprava do N//2)
+xf = fftfreq(N, dx)[:N//2]
+print("xf:", xf)
+# Amplituda (síla) signálu - musíme normalizovat dělením (2/N)
+amplitudy = (2.0/N) * np.abs(F_k[:N//2])
+print("Amplitudy:", amplitudy)
+# Ruční úprava nulté frekvence (DC složky - to je ta "1" na začátku rovnice)
+amplitudy[0] = amplitudy[0] / 2 
+
+# Vykreslíme to jako sloupečky (stem plot)
+plt.stem(xf, amplitudy)
+plt.xlim(-0.1, 2) # Přiblížíme začátek grafu, aby byly vidět špičky
+plt.title("Spektrum (FFT) - Všimni si špiček na frekvencích, které odpovídají rovnici")
+plt.xlabel("Frekvence (Hz)")
+plt.ylabel("Amplituda")
+plt.grid()
+
+plt.tight_layout()
+plt.show()

Old/1_fft_basics/1ext.py

@@ -0,0 +1,140 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.fft import fft, ifft, fftfreq
+p = False
+
+# ==========================================
+# 1. Definice mřížky (Grid) a signálu
+# ==========================================
+N = 100 
+x, dx = np.linspace(0, 2*np.pi, N, endpoint=False, retstep=True)
+
+# Rozložení signálu na jednotlivé komponenty (abychom viděli, co se z čeho skládá)
+f1 = 1 * np.ones_like(x)         # DC složka (konstanta 1)
+f2 = 2 * np.sin(x)               # Frekvence 1
+f3 = 10 * np.sin(5*x)            # Frekvence 5 (sinusovka)
+f4 = 3 * np.cos(5*x)             # Frekvence 5 (kosinusovka)
+
+# Výsledná složená funkce ze zadání
+f_x = f1 + f2 + f3 + f4
+
+if p:
+    print("dx:", dx)
+    print("x:", x)
+    print("f(x):", f_x)
+
+# ==========================================
+# 2. Výpočet FFT a zpětné transformace
+# ==========================================
+F_k = fft(f_x)
+f_x_reconstructed = ifft(F_k)
+
+# Výpočet frekvencí pro osu X
+xf_vsechny = fftfreq(N, dx)
+
+# Pro srozumitelné spektrum nás většinou zajímá jen kladná část spektra
+xf_kladne = xf_vsechny[:N//2]
+F_k_kladne = F_k[:N//2]
+
+# Amplitudy (síla signálu) - normalizace pomocí (2/N)
+amplitudy = (2.0/N) * np.abs(F_k_kladne)
+amplitudy[0] = amplitudy[0] / 2  # DC složka se nedělí dvěma, vracíme ji na správnou hodnotu
+
+# Fáze (posun signálu v radiánech)
+faze = np.angle(F_k_kladne)
+# Odfiltrování šumu z fáze: tam, kde je nulová nebo minimální amplituda, nemá smysl počítat fázi
+faze[amplitudy < 0.1] = 0
+
+# Výkonové spektrum (Power Spectrum - energie na dané frekvenci)
+vykon = np.abs(F_k_kladne)**2
+
+
+# ==========================================
+# 3. VYKRESLENÍ - OBRÁZEK 1: Časová oblast
+# ==========================================
+plt.figure(figsize=(12, 10))
+plt.subplot(3, 1, 1)
+plt.plot(x, f_x, label="Výsledný signál f(x)", color='black', linewidth=2)
+plt.title("1. Celkový složený signál f(x) v čase")
+plt.ylabel("Amplituda")
+plt.grid(True); plt.legend()
+
+plt.subplot(3, 1, 2)
+plt.plot(x, f1, label="DC (1)", linestyle='--')
+plt.plot(x, f2, label="2*sin(x)", linestyle='--')
+plt.plot(x, f3, label="10*sin(5x)", linestyle='--')
+plt.plot(x, f4, label="3*cos(5x)", linestyle='--')
+plt.title("2. Jednotlivé " + "skryté" + " složky, ze kterých se signál skládá")
+plt.ylabel("Amplituda")
+plt.grid(True); plt.legend()
+
+plt.subplot(3, 1, 3)
+plt.plot(x, f_x, label="Původní f(x)", linewidth=4, alpha=0.5)
+plt.plot(x, f_x_reconstructed.real, '--', label="Složená (iFFT)", color='red')
+plt.title("3. Původní signál vs. Rekonstruovaný signál (Úspěšný test!)")
+plt.xlabel("x (Čas / Prostor)")
+plt.ylabel("Amplituda")
+plt.grid(True); plt.legend()
+plt.tight_layout()
+
+
+# ==========================================
+# VYKRESLENÍ - OBRÁZEK 2: Surová komplexní čísla z FFT
+# ==========================================
+plt.figure(figsize=(12, 8))
+# Reálná část reprezentuje podíly KOSINUSOVÝCH funkcí
+plt.subplot(2, 1, 1)
+plt.stem(range(N), F_k.real, basefmt=" ")
+plt.title("4. Surový výstup FFT: Reálná část (Odpovídá zastoupení KOSINUSŮ ve frekvencích)")
+plt.ylabel("Reálná hodnota")
+plt.grid(True)
+
+# Imaginární část reprezentuje podíly SINUSOVÝCH funkcí
+plt.subplot(2, 1, 2)
+plt.stem(range(N), F_k.imag, basefmt=" ")
+plt.title("5. Surový výstup FFT: Imaginární část (Odpovídá zastoupení SINUSŮ ve frekvencích)")
+plt.xlabel("Index pole (frekvenční koš)")
+plt.ylabel("Imaginární hodnota")
+plt.grid(True)
+plt.tight_layout()
+
+
+# ==========================================
+# VYKRESLENÍ - OBRÁZEK 3: Polární forma (Amplituda a Fáze)
+# ==========================================
+plt.figure(figsize=(12, 8))
+plt.subplot(2, 1, 1)
+plt.stem(xf_kladne, amplitudy, basefmt=" ")
+plt.xlim(-0.1, 1.2) # Ořízneme graf jen na zajímavou část, ať vidíme špičky
+plt.title("6. Amplitudové spektrum - To nejdůležitější (Síla jednotlivých frekvencí)")
+plt.ylabel("Amplituda")
+plt.grid(True)
+
+# Přidání textu přímo nad špičky v grafu pro přehlednost
+for i, amp in enumerate(amplitudy):
+    if amp > 0.5:
+        plt.text(xf_kladne[i], amp + 0.5, f"{amp:.1f}", ha='center', fontweight='bold')
+
+plt.subplot(2, 1, 2)
+plt.stem(xf_kladne, faze, basefmt=" ")
+plt.xlim(-0.1, 1.2)
+plt.title("7. Fázové spektrum - Fázový posun odhalených vln (v radiánech)")
+plt.xlabel("Frekvence")
+plt.ylabel("Fáze [rad]")
+plt.grid(True)
+plt.tight_layout()
+
+
+# ==========================================
+# VYKRESLENÍ - OBRÁZEK 4: Výkonové spektrum (Power Spectrum)
+# ==========================================
+plt.figure(figsize=(12, 4))
+plt.stem(xf_kladne, vykon, basefmt=" ", markerfmt='ro', linefmt='r-')
+plt.xlim(-0.1, 1.2)
+plt.title("8. Výkonové spektrum (Power Spectral Density) - Kde je ukryta většina energie signálu")
+plt.xlabel("Frekvence")
+plt.ylabel("Výkon (Absolutní hodnota na druhou)")
+plt.grid(True)
+plt.tight_layout()
+
+plt.show()

Old/2_grf_1d/fft1d.py

@@ -0,0 +1,108 @@
+"""
+Generování náhodných polí pomocí spektrálních metod (FFT / NUFFT).
+Úkoly 3.3:
+  1) make_hermitian    – reálná funkce z náhodných komplexních koeficientů
+  2) generate_grf_fft  – GRF na pravidelné mřížce (Gaussovská korelace)
+  3) generate_grf_nufft – GRF na nepravidelné mřížce (finufft)
+
+    literatura
+
+"""
+
+import numpy as np
+from scipy.fft import fftfreq
+import finufft
+
+
+# SPEKTRÁLNÍ HUSTOTA
+
+def spectral_density_gaussian(omega, phi, sigma=1.0):
+    # C(r) = sigma^2 * exp(-r^2 / 2phi^2)  =>  S(w) = sigma^2 * sqrt(2pi)*phi * exp(-w^2*phi^2/2)
+    return sigma**2 * np.sqrt(2 * np.pi) * phi * np.exp(-0.5 * (omega * phi)**2)
+
+def spectral_density_exponential(omega, phi, sigma=1.0):
+    # C(r) = sigma^2 * exp(-|r|/phi)  =>  S(w) = sigma^2 * 2phi / (1 + (w*phi)^2)
+    return sigma**2 * 2 * phi / (1 + (omega * phi)**2)
+
+
+# HERMITOVSKÁ SYMETRIE – zaručí reálný výstup IFFT
+
+def make_hermitian(f_hat, centered=False):
+    """
+    Vyrobí hermitovsky symetrické koeficienty: f_hat[-k] = conj(f_hat[k])
+
+    centered=False  – FFT pořadí:      DC na indexu 0, záporné frekvence na konci
+    centered=True   – centrované pořadí: DC uprostřed na indexu N//2
+    """
+    N = len(f_hat)
+    f = f_hat.copy()
+
+    if centered:
+        zi = N // 2                           # index DC (nulové frekvence)
+        f[zi] = f[zi].real
+        for i in range(1, N // 2):
+            f[zi - i] = np.conj(f[zi + i])   # záporné frekvence = sdružené kladných
+    else:
+        f[0] = f[0].real                      # DC
+        if N % 2 == 0:
+            f[N // 2] = f[N // 2].real        # Nyquistova frekvence
+        for k in range(1, N // 2):
+            f[-k] = np.conj(f[k])
+
+    return f
+
+
+# FFT – PRAVIDELNÁ MŘÍŽKA
+
+def generate_grf_fft(L, N, phi, sigma=1.0, seed=None,
+                     spectral_density_fn=spectral_density_gaussian):
+    """
+    Generuje 1D náhodné pole na pravidelné mřížce.
+      - náhodné koeficienty z N(0,1) komplexního rozdělení
+      - modulace filtrem sqrt(S(omega))
+      - hermitovská symetrie -> reálný výstup přes IFFT
+    spectral_density_fn: spectral_density_gaussian | spectral_density_exponential
+    """
+    rng = np.random.default_rng(seed)
+    dx = L / N
+    x = np.arange(N) * dx
+
+    omega = 2 * np.pi * fftfreq(N, d=dx)  # FFT pořadí
+    S_k = spectral_density_fn(omega, phi, sigma)
+
+    white = (rng.standard_normal(N) + 1j * rng.standard_normal(N)) / np.sqrt(2)
+    f_hat = make_hermitian(white * np.sqrt(S_k), centered=False)
+
+    field = np.fft.ifft(f_hat).real * N / np.sqrt(L)
+    return x, field
+
+
+# NUFFT – NEPRAVIDELNÁ MŘÍŽKA
+
+def generate_grf_nufft(L, N_freq, M_points, phi, sigma=1.0, seed=None, x_points=None,
+                       spectral_density_fn=spectral_density_gaussian):
+    """
+    Generuje 1D náhodné pole na nepravidelné mřížce pomocí NUFFT typu 2
+    (pravidelné frekvence -> nepravidelné prostorové body).
+    finufft očekává souřadnice v [-pi, pi].
+    """
+    rng = np.random.default_rng(seed)
+
+    # centrované pořadí: k = -N/2, ..., -1, 0, +1, ..., N/2-1
+    k = np.arange(-N_freq // 2, N_freq // 2)
+    omega = k * (2 * np.pi / L)
+    S_omega = spectral_density_fn(omega, phi, sigma)
+
+    white = (rng.standard_normal(N_freq) + 1j * rng.standard_normal(N_freq)) / np.sqrt(2)
+    f_hat = make_hermitian(white * np.sqrt(S_omega), centered=True)
+
+    if x_points is None:
+        x_nufft = rng.uniform(-np.pi, np.pi, M_points)
+    else:
+        x_nufft = (np.asarray(x_points) / L) * 2 * np.pi - np.pi  # [0,L] -> [-pi,pi]
+
+    field = finufft.nufft1d2(x_nufft, f_hat.astype(np.complex128))
+    x_final = (x_nufft + np.pi) * (L / (2 * np.pi))               # [-pi,pi] -> [0,L]
+    return x_final, field.real
+
+

Old/2_grf_1d/main2.py

@@ -0,0 +1,35 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import fft1d as rf
+
+L, N, phi = 100.0, 1024, 3.0
+
+fig, axes = plt.subplots(2, 2, figsize=(14, 8))
+
+# --- Gaussovská ---
+x, f = rf.generate_grf_fft(L, N, phi, seed=42)
+axes[0][0].plot(x, f, lw=1)
+axes[0][0].set_title(f"FFT Gaussovská (φ={phi})")
+
+x_nu, f_nu = rf.generate_grf_nufft(L, 256, 500, phi, seed=42)
+idx = np.argsort(x_nu)
+axes[0][1].plot(x_nu[idx], f_nu[idx], '-', alpha=0.4, color='gray', lw=0.8)
+axes[0][1].scatter(x_nu, f_nu, s=12, c=f_nu, cmap='viridis', zorder=3)
+axes[0][1].set_title(f"NUFFT Gaussovská (φ={phi})")
+
+# --- Exponenciální ---
+x, f = rf.generate_grf_fft(L, N, phi, seed=42, spectral_density_fn=rf.spectral_density_exponential)
+axes[1][0].plot(x, f, lw=1)
+axes[1][0].set_title(f"FFT Exponenciální (φ={phi})")
+
+x_nu, f_nu = rf.generate_grf_nufft(L, 256, 500, phi, seed=42, spectral_density_fn=rf.spectral_density_exponential)
+idx = np.argsort(x_nu)
+axes[1][1].plot(x_nu[idx], f_nu[idx], '-', alpha=0.4, color='gray', lw=0.8)
+axes[1][1].scatter(x_nu, f_nu, s=12, c=f_nu, cmap='viridis', zorder=3)
+axes[1][1].set_title(f"NUFFT Exponenciální (φ={phi})")
+
+for ax in axes.flat:
+    ax.grid(True, alpha=0.3)
+
+plt.tight_layout()
+plt.show()