Estensimetro ed FFT

Una per applicare un filtro ai dati dell'estensimetro bassato su analisi FFT (vedi qui e precedenti)

Dati originali con curva di interpolazione polinomiale

 Detrend sottraendo la polinomiale

Analisi FFT

Il primo picco ha un periodo di circa 24 giorni. I successivi sono relativi al ciclo giornaliero ed armoniche di ordine superiore

applicando un filtro passa basso si ottiene il grafico filtrato invertendo la FFT

import pandas as pd

import numpy as np

from matplotlib import pyplot as plt

from sklearn.preprocessing import PolynomialFeatures

from scipy.fft import fft, ifft

# legge il csv e lo mette in un dataframe

dati = pd.read_csv('dati_small_prg.csv', sep=',', parse_dates=['Data'], index_col=['Data'])

# taglia la prima parte dei dati

dati = dati[1:8500]

# resample

#print("**********************************************")

#dati.info()

#dati.resample('60min').mean()

#dati.info()

#print("**********************************************")

####### regressione

import scipy

def f(x, a, b, c):

return a*x**2 + b*x + c

#return a*np.log(x)+b

#return a*x+b

def residual(p, x, y):

return y - f(x, *p)

p0 = [1., 1., 1.]

x = np.arange(1,8500,1)

popt, pcov = scipy.optimize.leastsq(residual, p0, args=(x,dati['Est.[mm]']))

print("Parametri della regressione ")

print(popt)

print("Errore sui parametri della regressione")

print(pcov)

yn = f(x, *popt)

plt.plot(x, yn, 'or')

plt.plot(x, dati['Est.[mm]'])

plt.show()

#Calcolo dell'errore quadratico medio

#tra i dati veri e la curva di correlazione

# come radice quadrata della somma delle differenze al quadrato

# diviso per il numero di elementi

differenza = yn-dati['Est.[mm]']

err = np.square(yn-dati['Est.[mm]'])

s_err = np.sqrt(np.sum(err, axis=0)/err.shape)

print("Errore quadratico medio (mm)")

print(s_err)

np.savetxt('correla.csv', yn, delimiter=',', fmt='%f')

# cancella le colonne che non servono

dati.drop(dati.columns[[0,1,4,5]], axis=1, inplace=True)

print(dati.info())

#salva i dati sul csv

dati.to_csv('./nuovo.csv')

#https://docs.scipy.org/doc/scipy/reference/tutorial/fft.html

#https://pythonnumericalmethods.berkeley.edu/notebooks/chapter24.04-FFT-in-Python.html

plt.plot(differenza, label='Detrend')

plt.legend(loc='best')

plt.show()

#con il passo di campionamento di 1 ora c'e' una frequenza

# di 0.0008333 Hz

# corrispondente ad una 1 misura ogni 20 minuti 1/(60x20)

dati_fft = np.fft.fft(differenza)

N = len(dati_fft)

n = np.arange(N)

# sr sampling rate

# n numero di misure

sr = 0.00083333

T = N/sr

print(T)

freq = n/T

plt.plot(freq[1:4250], np.abs(dati_fft[1:4250]))

plt.xlabel('Freq (Hz)')

plt.ylabel('FFT Amplitude |X(freq)|')

plt.show()

#filtra

sig_fft_filtered = dati_fft.copy()

cut_off = 0.00083

sig_fft_filtered[np.abs(freq) < cut_off] = 0

filtered = ifft(sig_fft_filtered)

plt.plot(filtered)

plt.xlabel('Tempo')

plt.ylabel('mm')

plt.show()

################################################

# 4.87927E-7 (23.72 giorni)

# 1.16571E-5 (23.82 ore)

# 2.32256E-5 (2 armonica del secondo dato)

# 3.48999E-5 (3 armonica del secondo dato)

# 4.6367E-5 (4 armonica del secondo dato)

# 5.8022E-5 (5 armonica del secondo dato)

# 6.9479E-5 (6 armonica del secondo dato)

# 8.1181E-5 (7 armonica del secondo dato)

# 9.2730E-5 (8 armonica del secondo dato)