Esponente di Hurst su torbidita' di acquiferi apuani

In questo post vengono analizzati i valori di torbidita' della sorgente del Frigido (periodo 13/02/2020 01/12/2025 dati orari ) e Carbonera (26/07/2006 21/06/2013 dati ogni 3 minuti ) legate all'acquiferi carsico delle Alpi Apuane

 

Torbidita' del tempo Frigido

 

Il parametro Esponente di Hurst e' stato introdotto Hurst, H. E. (1951). Long-term storage capacity of reservoirs. Transactions of the American Society of Civil Engineers, 116(1), 770-808. 

E' una misura dell'inerzia del sistema idrogeologico (ovvero capacita' di stoccaggio) e della complessita' della rete carsica. Per i valore di 0.82 e 0.85 trovati si ha che il sistema non e' casuale (H=0.5), non e' deterministicamente semplice ed il comportamento e' correlato su differenti scale di tempo ed e' un sistema che ha "memoria" del proprio passato

Non e' possibile prevedere in maniera completamente deterministica il valore esatto del parametro indagato  

Frigido

Con modelli deterministici si avra' un discreto successo su scale di tempo brevi ma non sulla previsione a lunga scadenza. Il passaggio da regime deterministico a frattale e' indicato come tempo di crossover e si individua dai punti di flesso sull'andamento dei punti. Per Carbonera (grafico sottostante) il punto di cross over e'  tra 1000 e 10000 unita' che tradotto nel passo di campionamento si arriva ad intervallo da le 5 e le 50 ore

Carbonera

Valori di H vicini ad 1 indicano persistenza nei valori (sia che siano massimi che minimi) e questo puo' incidere anche sul valore massimo e non solo sulle code della distribuzione

La memoria non si esaurisce mai ma ha una legge di decadimento come una potenze (in questo caso 2*0.84 -2 = -0.32 

 

con questa legge e l'esponente calcolate la correlazione ha un peso superiore al  10% fino a 1000 intervalli temporali (rho = 0.10^(1/-0.32)

Il fatto che il calcolo di H risulti sostanzialmente identico su serie con diverso passo di campionamento e' un indice della frattalita' del fenomeno osservato (e' indipendente dalla scala del tempo di osservazione)

 

il valore di Hurst prossimo ad 1 indica che e' probabile che la torbidita' aumenti nel prossimo lag (passo di campionamento) e' fortemente dipendente da cio' che accade nel lag precedente . Cio' sarebbe in contrasto con il concetto di massimo relativo della curva di scarica ma la statistica deve essere considerata valida al mantenersi delle condizioni al contorno e delle relative forzanti. Con un valore di 0.82 la coda di esaurimento sara' lunga

 

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

# --- 1. FUNZIONE DEFINIZIONE per l'ANALISI R/S MANUALE ---

def compute_rs_data(series):

    """

    Calcola la statistica Rango Riscalato (R/S) per una serie temporale.

    Restituisce gli intervalli temporali (n) e i valori medi di R/S.

    """

    N = len(series)

    # Definisce gli intervalli temporali (n) in progressione geometrica fino a N/4

    n_values = np.unique(np.logspace(1, np.log10(N // 4), num=50, dtype=int))

    n_values = n_values[n_values >= 1]

   

    rs_values = []

   

    for n in n_values:

        if n == 0: continue

       

        num_blocks = N // n

        rs_per_block = []

       

        for i in range(num_blocks):

            block = series[i * n : (i + 1) * n]

           

            mean = np.mean(block)

            std_dev = np.std(block)

           

            if std_dev == 0:

                continue

            X = np.cumsum(block - mean)

            R = np.max(X) - np.min(X)

           

            rs_per_block.append(R / std_dev)

       

        if rs_per_block:

            rs_values.append(np.mean(rs_per_block))

        else:

            break

    return n_values[:len(rs_values)], np.array(rs_values)

# --- 2. CARICAMENTO E PREPARAZIONE DEI DATI ADATTATI ---

FILE_NAME = "frigido.csv"

COLUMN_NAME = "Torbidita NTU"

FILTER_N_THRESHOLD = 100

try:

    # ATTENZIONE: Caricamento con separatore PUNTO E VIRGOLA (;) e decimale VIRGOLA (,)

    df = pd.read_csv(FILE_NAME, sep=';', decimal=',')

   

    # Seleziona la colonna 'Torbidita NTU'

    discharge_series = df[COLUMN_NAME]

   

    # Conversione in numerico e gestione dei NaN (valori mancanti o non validi)

    discharge_series = pd.to_numeric(discharge_series, errors='coerce').dropna()

   

    # Converti in array NumPy

    discharge_data = discharge_series.values

   

except FileNotFoundError:

    print(f"Errore: Il file '{FILE_NAME}' non è stato trovato. Assicurati che sia nella directory corretta.")

    exit()

except KeyError:

    print(f"Errore: La colonna '{COLUMN_NAME}' non è stata trovata nel file. Controlla il nome esatto.")

    print(f"Colonne disponibili: {df.columns.tolist()}")

    exit()

if len(discharge_data) < 100:

     print("Errore: Non ci sono abbastanza punti dati validi per eseguire l'analisi R/S.")

     exit()

# --- 3. CALCOLO R/S ---

rs_n, rs_rs = compute_rs_data(discharge_data)

# Calcolo di H (Esponente di Hurst) per l'INTERVALLO COMPLETO

log_n_full = np.log10(rs_n)

log_rs_full = np.log10(rs_rs)

coefficients_full = np.polyfit(log_n_full, log_rs_full, 1)

H_full = coefficients_full[0]

c_full = 10**coefficients_full[1]

# --- 4. PLOT: ANALISI R/S INTERVALLO COMPLETO ---

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

plt.plot(rs_n, c_full * rs_n**H_full, color="deepskyblue", label=f'Pendenza Frattale (H={H_full:.2f})')

plt.scatter(rs_n, rs_rs, color="navy", alpha=0.5, label='Dati R/S empirici')

plt.xscale('log')

plt.yscale('log')

plt.xlabel('Intervallo Temporale (n)')

plt.ylabel('Rango Riscalato R/S')

plt.title(f'Scaling Frattale della Torbidità (Full Range)')

plt.legend()

plt.savefig('Frigido_Full_R_S_Plot.png')

# --- 5. PLOT: INTERVALLO FILTRATO (n > 100) ---

filter_mask = rs_n > FILTER_N_THRESHOLD

X_filtered = rs_n[filter_mask]

Y_filtered = rs_rs[filter_mask]

# Calcolo di H per il solo intervallo filtrato

if len(X_filtered) > 1:

    log_n_filtered = np.log10(X_filtered)

    log_rs_filtered = np.log10(Y_filtered)

   

    coefficients_filtered = np.polyfit(log_n_filtered, log_rs_filtered, 1)

    H_filtered = coefficients_filtered[0]

    c_filtered = 10**coefficients_filtered[1]

   

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

    plt.scatter(X_filtered, Y_filtered, color="darkgreen", alpha=0.7, label=f'Dati R/S empirici (n > {FILTER_N_THRESHOLD})')

    plt.plot(X_filtered, c_filtered * X_filtered**H_filtered, color="red", linestyle='--',

             label=f'Pendenza Frattale Filtrata (H={H_filtered:.2f})')

    plt.xscale('log')

    plt.yscale('log')

    plt.xlabel('Intervallo Temporale (n)')

    plt.ylabel('Rango Riscalato R/S')

    plt.title(f'Scaling Frattale della Torbidità (Intervallo Temporale > {FILTER_N_THRESHOLD})')

    plt.legend()

    plt.savefig('Frigido_Filtered_R_S_Plot_gt_100.png')

   

    print("\n--- Risultati del Calcolo ---")

    print(f"Hurst Esponent (Torbidità - Intervallo Completo): {H_full:.4f}")

    print(f"Hurst Esponent (Torbidità - Filtrato n>{FILTER_N_THRESHOLD}): {H_filtered:.4f}")

else:

    print(f"Nota: Non ci sono abbastanza punti dati rimasti dopo il filtraggio (n > {FILTER_N_THRESHOLD}) per generare il grafico filtrato.")

 

 

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)

Connettore N (disavventure di un radioamatore)

Pensavo che le antenne per radioamatori avessero connettori SMA, BNC e SO 239 (connettore HF femmina) e PL259 (connettore HF maschio)

Comprando una antenna Diamond x50n ho scoperto che esistono i connettori N

 

 

 

Frattalita' in spessore di strato

Aggiormento: un tentativo piu' serio e' qui 

Sempre leggendo i libri dei post precedenti ho visto che e' stata trovata una correlazione frattale nello spessore degli strati turbiditici arenacei. Avendo fatto la tesi sulle Arenarie del Cervarola mi sono incuriosito ma ho misure di sezioni o affioramenti

Gli unici dati che ho trovato in merito a spessori di strato sono relativi agli stratotipo ed in particolare dell'affioramento di Punta Piccola per il Piacenziano (The Global Standard Stratotype-section and Point (GSSP) of the Piacenzian Stage (Middle Pliocene) by D. Castradori, D. Rio, F. J. Hilgen, and L. J. Lourens https://doi.org/10.18814/epiiugs/1998/v21i2/003.Non si tratta di una formazione turbiditica ma non avendo altri dati ci ho provato 

Questa e' la distribuzione di frequenza degli spessori..un limite e' che la distribuzione non e' oltre un ordine di grandezza 

 

 plottando log-log lo spessore degli strati con il numero di strati con spessore superiore alla soglia si ha che

 

usando lo spettro di potenza 

il valore β≈0.47 e' compreso tra il rumore bianco ed ed il fattore di scala 1/f dove β=1 (Rothman, D. H., Grotzinger, J. P., & Richter, F. M. (1994). Scaling in the Sedimentary Rock Record. Nature, 372(6505), 337–342. suggerisce che il valore atteso per le sequenze stratigrafiche e' β=1 altrimenti detto Pink Noise )

Cerchero' dati piu' accurati

 

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)