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)

Kiss FFT

------------------------------------------------------------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <kiss_fft.h>
#include <tools/kiss_fftr.h>


#define NUM_FFT 256

int main()
{

float samples[NUM_FFT];


// si crea un segnale fittizio
for (int i = 0; i < NUM_FFT; i++) {
    samples[i] = sin(i*0.1);
    // con questo incremento si che i varia da 0 a 25.6 rad
    // ovvero circa 4 volte 2pigreco (ovvero 4 circonferenze)


}

    int isinverse = 1;
    kiss_fft_scalar zero;
    memset(&zero,0,sizeof(zero));
    kiss_fft_cpx fft_in[NUM_FFT];
    kiss_fft_cpx fft_out[NUM_FFT];
    kiss_fft_cpx fft_reconstructed[NUM_FFT];
    kiss_fftr_cfg fft = kiss_fftr_alloc(NUM_FFT ,0 ,0,0);
    kiss_fftr_cfg ifft = kiss_fftr_alloc(NUM_FFT,isinverse,0,0);

// azzera le matrici
for (int i = 0; i < NUM_FFT; i++) {
    fft_in[i].r = zero;
    fft_in[i].i = zero;
    fft_out[i].r = zero;
    fft_out[i].i = zero;
    fft_reconstructed[i].r = zero;
    fft_reconstructed[i].i = zero;
}

//inserisce i campioni nella parte reale della matrice di input
for (int i = 0; i < NUM_FFT; i++) {
     fft_in[i].r = samples[i];
     fft_in[i].i = zero;
     fft_out[i].r = zero;
     fft_out[i].i = zero;
 }

 kiss_fftr(fft, (kiss_fft_scalar*) fft_in, fft_out);
 kiss_fftri(ifft, fft_out, (kiss_fft_scalar*)fft_reconstructed);

// calcolo della potenza
for (int i = 0; i < NUM_FFT/2; i++) {
    samples[i] = sqrt((fft_out[i].i*fft_out[i].i)+(fft_out[i].i*fft_out[i].i))/(NUM_FFT*2);
    printf("%.6f\n\r", samples[i]);

}

    printf("Terminato!\n\r");
    return 0;
}

------------------------------------------------------------------------------------------------

Dato input

Potenza

Segnale ricostruito

FFT su Arduino

Era da tanto che volevo provare ad elaborare i dati di un accelerometro in FFT con Arduino. Ora grazie a questo post ne sono venuto a capo.

Ho usato un accelerometro GY-61 con l'uscita Z connessa ad A0 ed uno schermo Oled I2C per visualizzare i dati (collegamento VCC 5V , GND, SDA -> A4, SCL -> A5).
Al momento di compilare il codice completo Arduino IDE ha segnalato che lo sketch era troppo grande ed ho dovuto eliminare l'output dei dati sulla porta seriale

Di seguito un video di esempio con la registrazione dei passi e battendo sul tavolo con la mano

--------------------------------------------------------------------------------------------
#include "arduinoFFT.h"

#include <SPI.h>
#include <Wire.h>
#include <Adafruit_GFX.h>
#include <Adafruit_SSD1306.h>

Adafruit_SSD1306 display(4);

#define SAMPLES 128             
#define AMPLIFICA 10

arduinoFFT FFT = arduinoFFT();

unsigned int sampling_period_us;
unsigned long microseconds;

double vReal[SAMPLES];
double vImag[SAMPLES];

void setup() {
    sampling_period_us = round(1000000*(1.0/100));
    display.begin(SSD1306_SWITCHCAPVCC, 0x3C);
}

void loop() {
   
    display.clearDisplay();
    
    for(int i=0; i<SAMPLES; i++)
    {
        microseconds = micros();    
        vReal[i] = analogRead(0);
        vImag[i] = 0;
        while(micros() < (microseconds + sampling_period_us)){
        }
    }

    /*FFT*/
    FFT.Windowing(vReal, SAMPLES, FFT_WIN_TYP_HAMMING, FFT_FORWARD);
    FFT.Compute(vReal, vImag, SAMPLES, FFT_FORWARD);
    FFT.ComplexToMagnitude(vReal, vImag, SAMPLES);

      
    for(int i=2; i<(SAMPLES/2)-2; i++)  //se si fa partire i da zero si prende anche la parte continua del segnale
    {
        display.drawLine(i, 0, i, AMPLIFICA*vReal[i]*(5.0/1023.0), WHITE);
    }
    display.display();
    delay(1);  
}

Motore acceso??

Il problema e' questo : e' possibile determinare se una macchina ferma ha il motore acceso o meno usando un telefono Android??

Si tratta di un problema per una app che dovevo fare ma che poi non ha mai visto la luce e la risposta e' si'....usando la FFT sui dati dell'accelerometro

Il primo grafico e' relativo ad una misura di 30 secondi a macchina ferma e spenta, il secondo sempre su una acquisizione di 30 secondi a macchina ferma ed accesa al minimo del motore

E' abbastanza evidente una frequenza di risonanza della macchina su tutti e tre gli assi intorno ai 40 Hz (il motore girava intorno ai 900 rpm ma quello che si registra e' la vibrazione di tutta la struttura)

FFT su Arduino

Usando sempre come input il microfono amplificato ho provato la libreria ArduinoFFT ed in particolare l'esempio FFT_03 modificando il numero di samples a 128 e la frequenza di campionamento a 9600


--------
/*

Example of use of the FFT libray to compute FFT for a signal sampled through the ADC.
        Copyright (C) 2017 Enrique Condes

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#include "arduinoFFT.h"

arduinoFFT FFT = arduinoFFT(); /* Create FFT object */
/*
These values can be changed in order to evaluate the functions
*/
#define CHANNEL A0
const uint16_t samples = 128; //This value MUST ALWAYS be a power of 2
double samplingFrequency = 9600;

unsigned int delayTime = 0;

/*
These are the input and output vectors
Input vectors receive computed results from FFT
*/
double vReal[samples];
double vImag[samples];

#define SCL_INDEX 0x00
#define SCL_TIME 0x01
#define SCL_FREQUENCY 0x02

void setup()
{
  if(samplingFrequency<=1000)
    delayTime = 1000/samplingFrequency;
  else
    delayTime = 1000000/samplingFrequency;
  Serial.begin(115200);
  Serial.println("Ready");
}

void loop()
{
  for(uint16_t i =0;i<samples;i++)
  {
    vReal[i] = double(analogRead(CHANNEL));
    if(samplingFrequency<=1000)
      delay(delayTime);
    else
      delayMicroseconds(delayTime);
  }
  /* Print the results of the sampling according to time */
  Serial.println("Data:");
  PrintVector(vReal, samples, SCL_TIME);
  FFT.Windowing(vReal, samples, FFT_WIN_TYP_HAMMING, FFT_FORWARD); /* Weigh data */
  Serial.println("Weighed data:");
  PrintVector(vReal, samples, SCL_TIME);
  FFT.Compute(vReal, vImag, samples, FFT_FORWARD); /* Compute FFT */
  Serial.println("Computed Real values:");
  PrintVector(vReal, samples, SCL_INDEX);
  Serial.println("Computed Imaginary values:");
  PrintVector(vImag, samples, SCL_INDEX);
  FFT.ComplexToMagnitude(vReal, vImag, samples); /* Compute magnitudes */
  Serial.println("Computed magnitudes:");
  PrintVector(vReal, (samples >> 1), SCL_FREQUENCY);
  double x = FFT.MajorPeak(vReal, samples, samplingFrequency);
  Serial.println(x, 6);
  while(1); /* Run Once */

}

void PrintVector(double *vData, uint8_t bufferSize, uint8_t scaleType)
{
  for (uint16_t i = 0; i < bufferSize; i++)
  {
    double abscissa;
    /* Print abscissa value */
    switch (scaleType)
    {
      case SCL_INDEX:
        abscissa = (i * 1.0);
break;
      case SCL_TIME:
        abscissa = ((i * 1.0) / samplingFrequency);
break;
      case SCL_FREQUENCY:
        abscissa = ((i * 1.0 * samplingFrequency) / samples);
break;
    }
    Serial.print(abscissa, 6);
    Serial.print(" ");
    Serial.print(vData[i], 4);
    Serial.println();
  }
  Serial.println();
}
--------

Valori piu' grandi di 128 nel numero di samples generano un errore di out_of_memory

Suonando su una tastiera musicale la quinta ottava ho acquisito i dati e trascitto il valore di massimo picco

Come si vede i valori di frequenza calcolati sono molto differenti da quelli teorici ma si osservi che il rapporto di frequenza tra due note consecutive distante due semitoni e' molto prossimo al valore teorico di 1.1225. Inoltre il rapporto tra la frequenza teorica e quella calcolata e' sempre di circa 2.1 quindi puo' trattarsi della prima armonica superiore

Eseguendo la stessa prova ma sulla 4° ottava i risultati sono i seguenti

Con questo sketch le classi di frequenza hanno una larghezza di 75 Hz...quindi non e' comunque un sistema da usare come accordatore

In generale e' comunque un sistema che funziona bene nell'ambito delle frequenze udibili

FFT in Android e Java/NetBeans

Per elaborare i dati derivanti dai sensori posti sul telefono puo' essere utile una elaborazione della trasformata in Fourier; per il caso in esame sara' impiegata la libreria JTransform

La libreria si inserisce normalmente all'interno del progetto di NetBeans od Android

Per l'esempio non sono stati usati dati reali ma quelli calcolati dalla funzione

10sin(x)+5sin(x/2(

Grafici da WolframAlpha

il codice Java (quello Android e' identico) e'
------------------------------

package javaapplication1;

import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D;


/**
 *
 * @author l.innocenti
 */
public class JavaApplication1 {
   

    /**
     * @param args the command line arguments
     */
    public static void main(String[] args) {
       int k;
        // numero dei dati 
        int N = 1002;

       DoubleFFT_1D dfft = new DoubleFFT_1D(N);

       double[] data = new double[N];
       double[] outputData = new double[N];

//crea il segnale
       for (k=0;k<N-1;k++) {
            data[k] = 10*Math.sin(Math.toRadians(k))+5*Math.sin(Math.toRadians(k/2));
            //System.out.println(k+";"+data[k]);
        }
     
       
        // si usa Real FFT perche' dati sono reali
        dfft.realForward(data);  
        
        // il risultato viene messo nella matrice dei risultati secondo il formato
        // data[0] = somma di tutti gli input
        // se il numero di dati e' pari a partire dall'indice data[2] si ha la parte reale nell'indice pari e 
        // la parte immaginaria nell'indice dispari
        // data[2*k] = Re[k], 0 <= k < n / 2
        // data[2*k+1] = Im[k], 0 < k < n / 2
        // altrimenti se i dati sono dispari
        // data[2*k] = Re[k], 0 <= k < (n+1)/2
        // data[2*k+1] = Im[k], 0 < k< (n-1)/2
        // per calcolare il valore di frequenza di ogni array 
        // Fs = valore di campionamento in Hz
        // N = numero di punti
        // k = indice della matrice
        // Frequenza = (k*Fs) / n
        //
        //La potenza dello spettro e' data da 
        //sqrt(parte_reale*parte_reale+parte_immaginaria*parte_immaginaria)
        
        if(data.length%2==0){
            for(int i = 0; i < data.length/2; i++){
                outputData[i]= (float) Math.sqrt((Math.pow(data[2*i],2))+(Math.pow(data[2*(i)+1], 2)));
            }
        }else{
            for(int i = 0; i < data.length/2-1; i++){
                outputData[i]= (float) Math.sqrt((Math.pow(data[2*i],2))+(Math.pow(data[2*i+1], 2)));
            }
        }
    }
}

-------------------------------------

Il risultato e' il seguente