HyperCoast Python Library

Hypercoast e' una libreria Python per leggere dati da formati iperspettrali e visualizzare spettri

L'uso migliore e' quello all'interno di Jupyter Lab

 

import hypercoast

filepath = "ang20210604t105418_rfl_v2z1_img"

ds = hypercoast.read_aviris(filepath)

print(ds)

m = hypercoast.Map()

print(m)

m

m.add_aviris(ds, wavelengths=[1000, 700, 400], vmin=0, vmax=0.2)

m.add("spectral")

display(m)

 

Aviris-NG Grosseto

 

 

 

 

Analisi MNF su spettri di riflettanza di plastica

Devo cerca di lavorare su spettri di riflettanza di plastica e la prima domanda e': quale sono le bande significative?

Sono partito dal dataset delle plastiche (36 in tutto) contenuti in USGS Spectral Library

Per estrarre le feature ho usato MNF Maximum Noise Fraction (al posto di PCA)  

Nel grafico sottostante sono plottate le prime 3 componenti MNF 

La terza componente (quella di solito in cui sono presenti i dati spettrali, la prima componente e' influenzata dalla ) e' stata utilizzata per addestrare un rete random tree per trovare quali siano le bande effettivamente significative

 Per finire un confronto tra gli spettri non elaborati e la selezione random forest 

 

 Le bande piu' significative estratte 

 ind           nm     peso

150          600     0.018671
1289        1739   0.018409
149          599     0.014763
143          593     0.013682
1998        2448   0.012601 

 

import numpy as np

import glob

import os

from scipy.signal import savgol_filter

from scipy.linalg import eigh

folder = "./plastica"

files = sorted(glob.glob(os.path.join(folder, "*.txt")))

arrays = []

for f in files:

    data = np.loadtxt(f, skiprows=1)

    if data.shape[0] > 2100:

        arrays.append(data[100:2100])

       

X = np.array(arrays)

X = np.where((X >= 0) & (X <= 1), X, 0)

print(f"Shape originale: {X.shape}")

# 1. Smoothing e calcolo rumore

X_smooth = savgol_filter(X, window_length=11, polyorder=2, axis=1)

noise_residue = X - X_smooth

# 2. Calcolo matrici di covarianza

noise_cov = np.cov(noise_residue, rowvar=False)

data_cov = np.cov(X, rowvar=False)

# Regolarizzazione (Ridge) alla diagonale della matrice di rumore.

# Questo garantisce che la matrice sia definita positiva.

eps = 1e-6 * np.trace(noise_cov) / noise_cov.shape[0]

noise_cov_reg = noise_cov + np.eye(noise_cov.shape[0]) * eps

# 3. Esecuzione MNF

eigenvalues, eigenvectors = eigh(data_cov, noise_cov_reg)

idx = np.argsort(eigenvalues)[::-1]

eigenvalues = eigenvalues[idx]

eigenvectors = eigenvectors[:, idx]

# 5. Proiezione

X_mnf = np.dot(X, eigenvectors)

print(f"Prime 5 componenti MNF (SNR più alto): {eigenvalues[:5]}")

import matplotlib.pyplot as plt

# L'indice 0 corrisponde alla prima colonna dopo l'ordinamento decrescente

first_mnf_loadings = eigenvectors[:, 0]

wavelengths = np.arange(350, 2500)

wl_tagliate = wavelengths[100:2100]

idx_max = np.argmax(np.abs(first_mnf_loadings))

print(f"La banda più importante è la numero {idx_max}, che corrisponde a {wl_tagliate[idx_max]} nm")

wls = np.arange(350, 2501)[100:2100]

fig, axs = plt.subplots(3, 1, figsize=(12, 15), sharex=True)

colors = ['#1f77b4', '#ff7f0e', '#2ca02c']

for i in range(3):

    # Estraiamo i loadings per la componente i-esima

    loadings = eigenvectors[:, i]

   

    axs[i].plot(wls, loadings, color=colors[i], lw=1.5)

    axs[i].set_title(f"Loadings MNF Componente {i+1} (Autovalore: {eigenvalues[i]:.2e})")

    axs[i].set_ylabel("Peso")

    axs[i].grid(True, alpha=0.3)

   

    # Evidenziamo l'area diagnostica dei carbonati (2300-2350 nm)

    axs[i].axvspan(2300, 2350, color='red', alpha=0.1, label='Zona Carbonati')

axs[2].set_xlabel("Lunghezza d'onda (nm)")

plt.tight_layout()

plt.show()

from sklearn.ensemble import RandomForestRegressor

import pandas as pd

y = X_mnf[:, 2]

rf = RandomForestRegressor(n_estimators=100, random_state=42)

rf.fit(X, y)

importances = rf.feature_importances_

wls = np.arange(350, 2501)[100:2100]

feature_results = pd.DataFrame({'Wavelength': wls, 'Importance': importances})

# Ordiniamo per importanza decrescente

top_features = feature_results.sort_values(by='Importance', ascending=False).head(10)

print("Le 10 lunghezze d'onda più influenti (Feature Importance):")

print(top_features)

# 5. Visualizzazione

plt.figure(figsize=(10, 5))

plt.plot(wls, importances, color='purple', label='Feature Importance')

plt.fill_between(wls, importances, color='purple', alpha=0.3)

plt.title("Rilevanza delle Bande Spettrali (Random Forest)")

plt.xlabel("Lunghezza d'onda (nm)")

plt.ylabel("Importanza Relativa")

plt.grid(True, alpha=0.3)

plt.show()

plt.figure(figsize=(12, 7))

plt.plot(wls, X.T, color='steelblue', alpha=0.3)

plt.title(f"Spettri di Riflettanza del Dataset ({X.shape[0]} campioni)")

plt.xlabel("Lunghezza d'onda (nm)")

plt.ylabel("Riflettanza (0-1)")

plt.grid(True, linestyle='--', alpha=0.5)

# Per evitare legende infinite, ne mettiamo una generica

plt.plot([], [], color='steelblue', label='Spettri campioni')

plt.legend()

plt.show()

# --- GRAFICO 1: SPETTRI DI RIFLETTANZA ---

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 12), sharex=True)

ax1.plot(wls, X.T, color='steelblue', alpha=0.2)

ax1.plot(wls, np.mean(X, axis=0), color='red', lw=2, label='Spettro Medio')

ax1.set_title("Spettri di Riflettanza (Dati Originali Puliti)")

ax1.set_ylabel("Riflettanza")

ax1.grid(True, alpha=0.3)

ax1.legend()

# --- GRAFICO 2: RILEVANZA DELLE BANDE (Feature Importance) ---

ax2.plot(wls, importances, color='darkgreen', lw=1.5)

ax2.fill_between(wls, importances, color='darkgreen', alpha=0.2)

ax2.set_title("Rilevanza delle Bande Spettrali (Random Forest Importance)")

ax2.set_ylabel("Importanza Relativa")

ax2.set_xlabel("Lunghezza d'onda (nm)")

ax2.grid(True, alpha=0.3)

plt.tight_layout()

plt.show()

Slideshow Mandelbrot

 Ho trovato a questo indirizzo un file csv in cui sono riportate le coordinate di finestre di zoom per l'insieme di Mandelbrot in zone interessanti (sono oltre 105.000 righe)

 

 Quale migliore occasione per creare uno slideshow di immagini dell'insieme di Mandelbrot

 

 

 

import csv

import numpy as np

from PIL import Image

import matplotlib.cm as cm

# Parametri

width, height = 640, 480

max_iter = 1000

colormap = cm.turbo  # puoi cambiare: viridis, inferno, magma...

def mandelbrot_vectorized(re_min, re_max, im_min, im_max, width, height, max_iter):

    # Griglia complessa

    real = np.linspace(re_min, re_max, width)

    imag = np.linspace(im_max, im_min, height)  # y invertito

    C = real[np.newaxis, :] + 1j * imag[:, np.newaxis]

    Z = np.zeros(C.shape, dtype=np.complex64)

    escape_iter = np.zeros(C.shape, dtype=np.float32)

    mask = np.ones(C.shape, dtype=bool)  # punti ancora in evoluzione

    for n in range(max_iter):

        Z[mask] = Z[mask]**2 + C[mask]

        escaped = np.abs(Z) > 2

        newly_escaped = escaped & mask

        if np.any(newly_escaped):

            # escape velocity continua

            escape_iter[newly_escaped] = n + 1 - np.log2(np.log2(np.abs(Z[newly_escaped])))

        mask &= ~escaped  # aggiorna mask

    # punti che non fuggono restano max_iter

    escape_iter[mask] = max_iter

    return escape_iter

# Leggi CSV e genera immagini

with open("windows.csv", newline='', encoding='utf-8') as csvfile:

    reader = csv.DictReader(csvfile)

   

    for idx, row in enumerate(reader):

        row = {k.strip(): v for k, v in row.items()}

        re_min = float(row['re_min'])

        re_max = float(row['re_max'])

        im_min = float(row['im_min'])

        im_max = float(row['im_max'])

        M = mandelbrot_vectorized(re_min, re_max, im_min, im_max, width, height, max_iter)

        # Normalizza escape velocity su 0-1

        M_norm = M / M.max()

        # Applica colormap

        M_color = (colormap(M_norm)[:, :, :3] * 1000).astype(np.uint8)

        interior_mask = (M == max_iter)

        M_color[interior_mask] = [0, 0, 0]

        img = Image.fromarray(M_color)

        img.save(f"mandelbrot_escape_vectorized_{idx+1}.png")

        print(f"Saved mandelbrot_escape_vectorized_{idx+1}.png")

 

 

 

 

 

Montagne frattali

 Negli anni 90 quando i miei amici matematici leggevano (Ri)creazioni al calcolatore di Dewdney (vedi anche pag. 107 del libro The Magic Machine a handbook of computer sorcery  A.K.Dewdney 1990 oppure Le Scienze n 222 1987 pag 98 e succ)  andava di moda programmare il proprio algoritmo di morphing e di  creazioni di paesaggi frattali

L'algoritmo che mi avevano spiegato era 

1) partire da un triangolo

2) cercare un punto casuale interno al triangolo 

3) uniire i vertici del triangolo con il punto casuale interno ed iterare sui nuovi tre triangoli

 All'epoca non ero riuscito a fare il mio codice...oggi ho deciso di chiedera all'AI. la cosa interessante del codice proposto e' l'algoritmo con cui viene cercato il punto casuale interno nella funzione random_point_inside. Vengono generati due numeri casuali compresi tra 0 e 1. Se la somma e' maggiore di 1 viene ridotta perche' il sistema esplora in realta' rettangoli, se la somma e' maggiore di 1 si e' selezionato un punto non nel triangolo

import numpy as np

import matplotlib.pyplot as plt

from matplotlib.tri import Triangulation

# ------------------------------------------------------

# 1. SCEGLI UN TRIANGOLO DI PARTENZA

# ------------------------------------------------------

# Triangolo grande in coordinate 2D

points = [

    np.array([0.0, 0.0]),

    np.array([1.0, 0.0]),

    np.array([0.5, 0.9])

]

triangles = [(0,1,2)]     # un singolo triangolo iniziale

points = list(points)

# ------------------------------------------------------

# 2. FUNZIONE: scegli un punto interno casuale (barycentric random)

# ------------------------------------------------------

def random_point_inside_triangle(a,b,c):

    r1 = np.random.rand()

    r2 = np.random.rand()

    # riflessione per distribuzione uniforme

    if r1 + r2 > 1:

        r1 = 1 - r1

        r2 = 1 - r2

    return a + r1*(b-a) + r2*(c-a)

# ------------------------------------------------------

# 3. ITERAZIONE DI SUDDIVISIONE

# ------------------------------------------------------

def subdivide_random(triangles, points, iterations=500):

    for _ in range(iterations):

        # scegli casualmente un triangolo da suddividere

        idx = np.random.randint(len(triangles))

        i,j,k = triangles[idx]

        a,b,c = points[i], points[j], points[k]

        # nuovo punto interno

        p = random_point_inside_triangle(a,b,c)

        p_index = len(points)

        points.append(p)

        # crea tre nuovi triangoli

        new_tris = [

            (i, j, p_index),

            (j, k, p_index),

            (k, i, p_index)

        ]

        # sostituisci il triangolo originale con i 3 nuovi

        triangles[idx] = new_tris[0]

        triangles.extend(new_tris[1:])

    return triangles, points

triangles, points = subdivide_random(triangles, points, iterations=1500)

# ------------------------------------------------------

# 4. PLOT FINALE

# ------------------------------------------------------

pts = np.array(points)

tri = np.array(triangles)

plt.figure(figsize=(7,7))

plt.tripcolor(pts[:,0], pts[:,1], tri, edgecolors='k', linewidth=0.1, alpha=0.9)

plt.title("Random Triangle Subdivision – Paesaggio frattale 2D")

plt.axis("equal")

plt.show()

 

Un altro algoritmo (diamond square) e' basato su una griglia regolare lavorando sui vertici. Il modello e' frattale peri il metodo con cui viene calcolata la altezza  

 

 

 

 

 import numpy as np
import matplotlib.pyplot as plt

def midpoint_displacement(size, roughness=0.7):
    """
    Genera un terreno frattale 2D (matrice di elevazione Z)
    utilizzando l'algoritmo Diamond-Square semplificato (Midpoint Displacement).
    
    :param size: Il lato della griglia, deve essere (2^n) + 1.
    :param roughness: Controlla la "frattalità" (0=liscio, alto=irregolare).
    :return: Una matrice NumPy di elevazioni (Z).
    """
    # Inizializza la griglia Z con zeri
    Z = np.zeros((size, size), dtype=float)
    
    # 1. Inizializzazione degli angoli
    # Assegna valori casuali iniziali ai quattro angoli per dare il via al processo
    Z[0, 0] = np.random.rand() * 10
    Z[0, size - 1] = np.random.rand() * 10
    Z[size - 1, 0] = np.random.rand() * 10
    Z[size - 1, size - 1] = np.random.rand() * 10
    
    step = size - 1
    rand_range = 10.0  # Range iniziale della casualità
    
    # Processo iterativo
    while step > 1:
        half = step // 2
        
        # Quadrati (Diamond Step)
        # Calcola il punto medio di ogni quadrato (i, j)
        for i in range(0, size - 1, step):
            for j in range(0, size - 1, step):
                
                # Media dei quattro angoli del quadrato
                avg = (Z[i, j] + Z[i + step, j] + Z[i, j + step] + Z[i + step, j + step]) / 4.0
                
                # Aggiunge casualità al centro del quadrato
                Z[i + half, j + half] = avg + (np.random.rand() - 0.5) * 2 * rand_range
        
        # Diamanti (Square Step)
        # Calcola i punti medi dei bordi (punti 'diamante')
        for i in range(0, size - 1, half):
            for j in range((i + half) % step, size - 1, step):
                
                # Punti da mediare (gestione bordi per non sforare l'array)
                punti = [
                    Z[(i - half + size) % size, j], # Sopra
                    Z[(i + half) % size, j],         # Sotto
                    Z[i, (j - half + size) % size], # Sinistra
                    Z[i, (j + half) % size]          # Destra
                ]
                
                # Media dei quattro vicini (quelli che esistono)
                avg = np.mean(punti)
                
                # Aggiunge casualità al punto medio del bordo
                Z[i, j] = avg + (np.random.rand() - 0.5) * 2 * rand_range
                
                # Questo codice semplificato è per i bordi interni; un'implementazione
                # completa dovrebbe gestire meglio gli indici di confine
        
        # Riduce l'ampiezza di casualità per la prossima iterazione (l'algoritmo frattale)
        rand_range *= 2**(-roughness)
        step = half
        
    return Z

def visualize_wireframe(Z):
    """
    Visualizza la matrice di elevazione Z come un wireframe 3D.
    """
    size = Z.shape[0]
    
    # Creazione delle coordinate X e Y
    x = np.linspace(0, 1, size)
    y = np.linspace(0, 1, size)
    X, Y = np.meshgrid(x, y)

    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')

    # Plot del wireframe
    ax.plot_wireframe(X, Y, Z, color='brown', linewidth=0.3)

    # Impostazioni estetiche
    ax.set_title(f'Montagna Frattale Wireframe (Roughness={ROUGHNESS})')
    ax.set_xlabel('Asse X')
    ax.set_ylabel('Asse Y')
    ax.set_zlabel('Elevazione Z')
    
    # Regola gli assi Z per centrare meglio la montagna
    z_min, z_max = np.min(Z), np.max(Z)
    z_range = z_max - z_min
    ax.set_zlim(z_min - z_range * 0.1, z_max + z_range * 0.1)

    plt.show()

# --- Parametri di Esecuzione ---
# La dimensione (size) DEVE essere (2^n) + 1, ad esempio: 33, 65, 129, 257.
# 65 è un buon compromesso tra dettaglio e velocità.
SIZE = 65 
ROUGHNESS = 0.85 # Valori più alti (es. 0.9) danno montagne più "spigolose"

# 2. Generazione del Frattale e Visualizzazione
elevation_data = midpoint_displacement(SIZE, ROUGHNESS)
visualize_wireframe(elevation_data)

Frattalita' in misure geofisiche di pozzo (2)

Nella conclusione del post precedente era che i vari tentativi di ricondurre i dati sperimentali all'analisi frattale davano risultati contradittori...tutto cio' fino sono arrivato a leggere Fractal Behaviour of the Earth System V.P.Dimitri  Springer 2005 pag 63-81 Regularity Analysis Applied to Well Log Data (M.Fedi, D.Fiore, Mauro La Mann) dove parlando proprio dei dati di KST viene introdotto il concetto di multifrattalita' (per essere precisi di multifrattalita' aveva gia' parlato Mandelbrot ed in D. Turcotte  Fractals and Chaos in Geology and Geophysics si trova applicazione dei multifrattali applicata alle geoscienze a pag)

Prima di tutto si puo' effettuare un test se i dati del pozzo possono essere inquadrati all'interno della multifrattalita' 

Questo puo' essere fatto l'analisi di regolarita' locale tramite analisi di fluttuazione a momenti multipli (Multifractal Detrended Fluctuation Analysis, MF-DFA) 

Vengono testati valori del parametro q da -5 a 5 (escludendo lo 0) verificando la pendenza h(q) 

 

 

Se il segnale fosse monofrattale la curva sarebbe una linea orizzontale ma non e' il caso. Quindi il test mostra la presenza di trend multifrattale con q>0 (destra) influenzato dalle grandi fluttuazioni (picchi del grafico dei valori di Gamma o contrasti gelogici) o q<0 (sinistra) influenzati da regioni piu' omogenee o da rumore di fondo

La linea rossa (q=2) indica  la persistenza media dell'intera serie

nel grafico sottostante viene descritta l'intensita' della multifrattalita' (se non fossero dati multi frattali si ridurrebbe ad un punto. Piu' ampia e' l'apertura della curva maggiore e' l'intensita' della frattalita'

I dati quindi non possono essere descritti con un singolo esponente di Hurst 

 

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from scipy.stats import linregress

# ======================================================

# 0 — FUNZIONE MF-DFA MANUALE

# ======================================================

def mfdfa(series, qs, scales, order=1):

"""

Calcola l'esponente di Hurst generalizzato h(q) utilizzando l'algoritmo MF-DFA.

"""

N = len(series)

# 1. Integrazione del segnale (necessaria per il detrending)

Y = np.cumsum(series - np.mean(series))

h_q = []

for q in qs:

F_q_2 = [] # Momenti di fluttuazione di q-esimo ordine

for s in scales:

s = int(s)

# Suddivisione in segmenti di lunghezza s

segments = N // s

# Calcolo della varianza per tutti i segmenti

F_q_s_list = []

for i in range(segments):

idx = slice(i * s, (i + 1) * s)

Y_segment = Y[idx]

# 2. Detrending locale (regressione polinomiale di ordine 'order')

x_local = np.arange(s)

p = np.polyfit(x_local, Y_segment, order)

Y_fit = np.polyval(p, x_local)

# 3. Calcolo della varianza locale (residui al quadrato)

variance = np.mean((Y_segment - Y_fit)**2)

F_q_s_list.append(variance)

# 4. Calcolo della funzione di fluttuazione di q-esimo ordine

F_q_s_list = np.array(F_q_s_list)

if q != 2.0:

# Esecuzione del momento q-esimo (q/2 sulle varianze)

F_q = np.mean(F_q_s_list**(q / 2.0))**(1.0 / q)

else:

# Caso speciale q=2 (momento quadratico)

F_q = np.sqrt(np.mean(F_q_s_list))

F_q_2.append(F_q)

# 5. Regressione log-log: F_q(s) ~ s^h(q)

log_s = np.log10(scales)

log_F_q = np.log10(F_q_2)

# Filtra i valori infiniti o NaN che possono emergere da log(0)

valid_indices = np.isfinite(log_s) & np.isfinite(log_F_q)

if np.sum(valid_indices) > 1:

slope, _, _, _, _ = linregress(log_s[valid_indices], log_F_q[valid_indices])

h_q.append(slope)

else:

h_q.append(np.nan)

return np.array(h_q)

# ======================================================

# 1 — LETTURA FILE CSV E PREPARAZIONE DATI

# ======================================================

filename = "well_log.csv"

RIGHE_LIMITE = 12700

COLONNA_INDICE = 1

try:

df = pd.read_csv(filename, nrows=RIGHE_LIMITE)

except FileNotFoundError:

print(f"Errore: File '{filename}' non trovato.")

exit()

x = df.iloc[:, COLONNA_INDICE].values

x = pd.to_numeric(x, errors='coerce')

x = x[~np.isnan(x)]

# Riscalamento (opzionale)

xmin = np.min(x)

if xmin <= 0:

shift = abs(xmin) + 1e-6

x = x + shift

print(f"Riscalamento effettuato: shift = {shift:.6f}")

N = len(x)

print(f"Analisi eseguita su {N} campioni (Prime {RIGHE_LIMITE} righe del file).")

# ======================================================

# 2 — PARAMETRI MF-DFA

# ======================================================

# Scelta dell'intervallo q

Q_MIN = -5

Q_MAX = 5

Q_STEP = 0.5

qs = np.arange(Q_MIN, Q_MAX + Q_STEP, Q_STEP)

qs = qs[qs != 0] # Esclude q=0

# Scegliere un intervallo di scale logaritmico per una migliore regressione

min_scale = int(N / 50) # Scale minime per un detrending efficace (p.es. 50 campioni)

max_scale = int(N / 4) # Scala massima (1/4 della lunghezza totale)

scales = np.unique(np.logspace(np.log10(min_scale), np.log10(max_scale), num=50).astype(int))

ORDER = 1 # Ordine del polinomio di detrending locale

# ======================================================

# 3 — ESECUZIONE MF-DFA E CALCOLO h(q)

# ======================================================

print("Calcolo degli esponenti h(q)...")

hq = mfdfa(x, qs, scales, order=ORDER)

# Filtra i valori NaN (dovuti a una regressione fallita)

valid_hq = hq[np.isfinite(hq)]

valid_qs = qs[np.isfinite(hq)]

# ======================================================

# 4 — ANALISI E INTERPRETAZIONE DEI RISULTATI

# ======================================================

# 4a. Analisi monofattale (q=2)

if 2.0 in valid_qs:

q_mono_index = np.where(valid_qs == 2.0)[0][0]

hurst_monofractal = valid_hq[q_mono_index]

else:

# Interpolazione per q=2 se non è stato calcolato

hurst_monofractal = np.interp(2.0, valid_qs, valid_hq)

# 4b. Analisi della Multifrattalità

h_max = np.max(valid_hq)

h_min = np.min(valid_hq)

Delta_h = h_max - h_min

# ======================================================

# 5 — PLOT DEI RISULTATI

# ======================================================

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=False)

# --- A) Grafico h(q) vs q (Test di Multifattalità)

ax1.plot(valid_qs, valid_hq, 'o-', color='darkblue', markersize=4)

ax1.axhline(hurst_monofractal, color='red', linestyle='--',

label=f'H monofattale (q=2): {hurst_monofractal:.3f}')

ax1.set_title('Esponente di Hurst Generalizzato $h(q)$')

ax1.set_xlabel('Momento di Fluttuazione $q$')

ax1.set_ylabel('Esponente di Scaling $h(q)$')

ax1.grid(True)

ax1.legend()

# --- B) Lo Spettro Multifrattale f(alpha) vs alpha (Singolarità)

# Calcolo dello Spettro tramite Trasformata di Legendre

tau_q = valid_qs * valid_hq - 1

alpha = np.gradient(tau_q, valid_qs)

f_alpha = valid_qs * alpha - tau_q

ax2.plot(alpha, f_alpha, 's-', color='green', markersize=4)

ax2.set_title('Spettro Multifrattale $f(\\alpha)$')

ax2.set_xlabel('Esponente di Singolarità $\\alpha$')

ax2.set_ylabel('Dimensione Frattale $f(\\alpha)$')

ax2.grid(True)

plt.tight_layout()

plt.show()

print("\n" + "=" * 50)

print(f"ANALISI MULTIFRATTALE (MF-DFA) - RISULTATI")

print(f"Δh = h_max - h_min = {Delta_h:.4f}")

if Delta_h > 0.1:

print("CONCLUSIONE: Il segnale è **fortemente multifrattale** (Δh > 0.1).")

print("La sua complessità non è uniforme: sono presenti strutture a diverse scale e diversi tipi di singolarità (anomalie e omogeneità).")

else:

print("CONCLUSIONE: Il segnale è prevalentemente **monofattale** (Δh <= 0.1).")

print(f"Esponente di Hurst Monofattale (q=2): H ≈ {hurst_monofractal:.3f}")

print("=" * 50)