readerCoerente 0.9.0

Package to load files saved by Cohaerentia Interferometric Interrogator.

Sensore Coerente Cohaerentia - lettore file TDMS

Pacchetto per lettura e visualizzazione delle misure salvate da software per sensore Coerente sviluppato da Cohaerentia.

Installazione

Il pacchetto è disponibile su PyPi e può essere installato tramite pip:

pip install readerCoerente

In alternativa, scaricare il wheel dal repository e installare tramite pip:

pip install <path_to_wheel>

Utilizzo

Importare pacchetto:

from readerCoerente import File, multipleFiles, detrend,

Caricamento dei dati

Singolo file

Si utilizza la classe File, specificando il percorso al file da caricare.

prova = File('data/2023-07-19_15-15-20_prova0_Ph.tdms')
print(prova)

Dove si è usato print per ottenere alcune informazioni di base sul file caricato.

Multipli files concatenati

E' possibile concatenare i dati contenuti in multipli files, se contenuti nella stessa cartella folder e con le stesse proprietà (per esempio, stessa frequenza di campionamento). Il pacchetto non verifica la continuità temporale dei dati contenuti nei files. Indicare nel parametro filename il path "comune" nei files. Per esempio, se i files sono:

data/2023-07-19_15-23-53_prova3_Ph.tdms
data/2023-07-19_15-24-06_prova3_Ph.tdms
data/2023-07-19_15-24-12_prova3_Ph.tdms
data/2023-07-19_15-24-37_prova3_Ph.tdms

è sufficiente indicare filename = 'prova3'.

provaMulti = multipleFiles(folder="data",filename="prova3")
print(provaMulti)

E' possibile selezionare un intervallo temporale specifico per l'analisi, indicando i parametri start_time_utc e stop_time_utc (in formato datetime). Il default di start_time_utc è il primo istante di campionamento del primo file, mentre il default di stop_time_utc è l'ultimo istante di campionamento dell'ultimo file.

from datetime import datetime, timedelta
start_time = datetime(2025,9,8,10,0,0)
stop_time = start_time + timedelta(hours=6)

folder = "<folder_path>"
prova_range = multipleFiles(folder=folder,filename="<file_common_name>.tdms",
                            start_time_utc=start_time,stop_time_utc=stop_time)

Un esempio di caricamento di dati a blocchi, della durata di una settimana, partendo dal primo lunedì disponibile:

from glob import glob
from os.path import join
from datetime import datetime, timedelta
import numpy as np

def next_weekday(date: np.datetime64, weekday: int) -> np.datetime64:
    days_ahead = weekday - date.tolist().date().weekday()
    if days_ahead <= 0: # Target day already happened this week
        days_ahead += 7
    next_weekday_datetime = datetime.combine(date.tolist().date() + timedelta(days_ahead), datetime.min.time())
    return np.datetime64(next_weekday_datetime)

folder = "<folder_path>"
files = sorted(glob(join(folder, "<file_common_name>.tdms")))

first_time_available = File(files[0]).start_time
print("First time available:", first_time_available)
first_actual_monday = next_weekday(first_time_available, 0)
print("First actual monday:", first_actual_monday)
temp = File(files[-1])
last_time_available = temp.start_time + np.timedelta64(int(temp.n_samples / temp.freq_sampling * 1000), 'ms')
print("Last time available:", last_time_available)
measurement_duration = last_time_available - first_actual_monday
print("Measurement duration (days):", measurement_duration.astype('timedelta64[D]'))

days_per_plot = 7
number_of_plots = int(np.ceil(measurement_duration.astype('timedelta64[D]') / np.timedelta64(days_per_plot, 'D')))
print("Number of plots:", number_of_plots)

time_windows: list[tuple[np.datetime64, np.datetime64]] = list()
for plot_index in range(number_of_plots):
    start_time = first_actual_monday + np.timedelta64(plot_index * days_per_plot, 'D')
    end_time = min(start_time + np.timedelta64(days_per_plot, 'D'), last_time_available)
    print(f"Plot {plot_index+1}/{number_of_plots}: from {start_time} to {end_time}")
    dati = multipleFiles(folder=folder, filename="installazione", start_time_utc=start_time, stop_time_utc=end_time)

Lettura file .ini con parametri di normalizzazione

from readerCoerente import read_config
para = read_config("Cohaerentia3PD_003.ini")
print(para)

Metodi

Per accedere al numpy array della fase, per esempio per plottarlo:

import plotly.graph_objects as go
go.Figure(layout=dict(title=prova.filename, xaxis_title='Time (s)', yaxis_title='Phase (rad)'),
          data=[go.Scatter(x=prova.time, y=prova.phase)]).show()

Per convertire i radianti ad allungamento relativo misurato si può utilizzare il metodo convertToElongation che converte i radianti in metri. Di default viene utilizzato il coefficiente teorico per fibra in silica: rad = 9.239e6 metri. E' possibile utilizzarne uno differente attraverso il parametro opzionale coefficient.

allungamento = prova.convert_to_elongation()

O, analogamente, per convertire a variazione di temperatura, utilizzare il metodo .convert_to_temperature() che utilizza il coefficiente di temperatura per fibra pari a 42.56 rad/°C*m ed è quindi necessario indicare anche la lunghezza del ramo dell'interferometro. Anche in questo caso è possibile utilizzare un coefficiente differente specificandolo nel parametro coefficient.

Se presenti nel file TDMS, è possibile accedere a tutti i canali, per esempio per plottarne una parte:

for ch in prova.channels_name:
    data = getattr(prova,ch)

Analisi spettrali

Per calcolare la power spectrum density, per esempio della fase:

from scipy import signal

freqs, psd = signal.welch(sig, fs=prova.freq_sampling,nperseg=1024)

go.Figure(layout=dict(title='PSD: power spectral density',
                      xaxis_title='Frequency [Hz]',
                      yaxis_title='Power', yaxis_type='log'),
          data=[go.Scatter(x=freqs, y=psd)]).show()

Per calcolare la Short-Time Fourier Transform (STFT), si può utilizzare il metodo .stft() della classe File o multipleFiles.

frequencies,times,power = prova.stft( window_size=1024, dB_scale=True)

from plotly.subplots import make_subplots
import numpy as np

fig_stft = make_subplots(rows=2, cols=1, shared_xaxes=True,
                            row_heights=[0.3, 0.7],vertical_spacing=0.05)

fig_stft.add_trace(go.Heatmap(
    x=times,
    y=frequencies,
    z=power,  # Convert to dB scale (avoid log(0))
    colorbar=dict(title='Power (dB)'),
    colorscale='Jet',
),row=2, col=1)

fig_stft.add_trace(go.Scatter(
    x=prova.time,
    y=prova.phase,
    line=dict(color='blue')
), row=1, col=1)

fig_stft.update_layout(
    title=prova.filename,
    showlegend=False,
    yaxis_title='Phase (rad)',
    yaxis2_title='Frequency (Hz)',
    xaxis2_title='Time (s)',
    width=1000,
    height=800
)

Detrend

Per esempio, utilizzando la funzione detrend del pacchetto, che è un filtfilt passa alto.

ph_detrend = detrend(prova.phase, prova.freq_sampling, highpass_cutoff_Hz=100)

go.Figure(layout=dict(title=f"{prova.filename} detrended", xaxis_title='Time (s)', yaxis_title='Phase (rad)'),
          data=[go.Scatter(x=prova.time, y=ph_detrend)]).show()

Esportazione dei dati in file h5

Per visualizzare i file h5 si può utilizzare Visual Studio Code con l'estensione H5Web.

prova.export_to_h5('data/prova.h5')

Sottocampionamento

Il metodo della classe si chiama .undersampling() e accetta parametro opzionale chiamato factor il cui valore di default è pari a 100. Esiste poi parametro opzionale what che accetta valori {'all', 'phase', 'channel'} e serve per stabilire cosa viene sottocampionato.

Nel caso in cui si sia interessati a non perdere per sempre quei campioni sarebbe necessario o ricaricare il file, creando di fatto un nuovo oggetto da zero, oppure si deve copiare in un nuovo oggetto prima di effettuare il sottocampionamento. La copia di oggetti non può essere effettuata con un semplice:

provaunder = prova

ma è necessario utilizzare il pacchetto copy:

import copy
provaunder = copy.deepcopy(prova)

Asse dei tempi assoluta

Attenzione che plotly diventa particolarmente lento in questo caso; si consiglia quindi di sottocampionare i dati prima di plottarli. Utilizzando la proprietà .absolute_time:

import copy
under_data = copy.deepcopy(misure[0])
under_data.undersampling(factor=50)
go.Figure(data=go.Scatter(x=under_data.absolute_time, y=under_data.phase)).show()

Compatibilità con Obspy

Il pacchetto readerCoerente è compatibile con Obspy, in particolare con la classe Trace. E' possibile convertire un oggetto File o multipleFiles in un oggetto Trace.

from obspy.core.trace import Trace,Stats

dati = multipleFiles(folder=folder, filename="installazione")
print(f"Loaded data from {dati.start_time} to {dati.start_time + timedelta(seconds=len(dati.phase)/dati.freq_sampling)}")

obs_stats = Stats(header=dict(delta=1/dati.freq_sampling, starttime=dati.start_time.astype(datetime)))
obs_trace = Trace(data=dati.phase, header=obs_stats)
obs_trace.filter("highpass", freq=2,)

fig_dayplot = obs_trace.plot(#starttime=UTCDateTime("2025-09-28T18:00:00"),
    type='dayplot',interval=60,
    tick_format='%d/%m %H:%M', right_vertical_labels=True,
    title=f"{create_window_time_str(time_window[0], time_window[1])} - Detrended Phase (rad)", size=(2500,2000),
    outfile='plot.png', dpi=200)