mercoledì 11 febbraio 2026

Mappa NDVI medio annuale con Google Earth Engine

Al contrario del precedente post https://debiaonoldcomputers.blogspot.com/2026/02/analisi-multitemporale-ndvi-con-google.html qui l'idea e' di creare una mappa di NDVI al posto della serie tempo per vedere anomalie spaziali

 


 

Landsat 8 

 // 1. Define the Study Area

var areaOfInterest = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[10.904865172222156, 43.01615135702198],
          [10.904865172222156, 43.00648615894694],
          [10.925807860210437, 43.00648615894694],
          [10.925807860210437, 43.01615135702198]]], null, false);

// 2. Pre-processing Function for Landsat 8
var processL8 = function(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var ndvi = opticalBands.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');
  return image.addBands(ndvi).select('NDVI');
};

// 3. Load 2025 Median NDVI
var ndvi2025 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
  .filterBounds(areaOfInterest)
  .filterDate('2025-01-01', '2025-12-31')
  .map(processL8)
  .median()
  .clip(areaOfInterest);

// 4. Visualization Parameters
var palette = [
  '#FFFFFF', '#CE7E45', '#FCD163', '#66A000', '#207401', '#056201', '#011301'
];
var ndviViz = {min: 0, max: 0.9, palette: palette};

// 5. Display the Map
Map.addLayer(ndvi2025, ndviViz, 'NDVI Median 2025');

// --- 6. ADDING THE LEGEND ---

// Create the legend panel
var legend = ui.Panel({
  style: {
    position: 'bottom-left',
    padding: '8px 15px'
  }
});

// Create legend title
var legendTitle = ui.Label({
  value: 'NDVI Index (2025)',
  style: {fontWeight: 'bold', fontSize: '16px', margin: '0 0 4px 0', padding: '0'}
});
legend.add(legendTitle);

// Helper function to create rows in the legend
var makeRow = function(color, name) {
  var colorBox = ui.Label({
    style: {
      backgroundColor: color,
      padding: '8px',
      margin: '0 0 4px 0'
    }
  });
  var description = ui.Label({
    value: name,
    style: {margin: '0 0 4px 6px'}
  });
  return ui.Panel({
    widgets: [colorBox, description],
    layout: ui.Panel.Layout.Flow('horizontal')
  });
};

// Define labels for the palette colors
var names = ['Low/No Veg (< 0.1)', 'Soil / Sparse', 'Moderate Veg', 'Dense Veg', 'Healthy Forest', 'Peak Vigour (> 0.8)'];
var colors = ['#FFFFFF', '#CE7E45', '#FCD163', '#66A000', '#207401', '#011301'];

// Add rows to legend
for (var i = 0; i < 6; i++) {
  legend.add(makeRow(colors[i], names[i]));
}

// Add legend to the map
Map.add(legend);

Export.image.toDrive({
  image: ndvi2025,            // The variable name of your NDVI map
  description: 'NDVI_2025_Tuscany', 
  folder: 'EarthEngine_Exports', // Name of the folder in your Google Drive
  fileNamePrefix: 'NDVI_2025_Site',
  region: areaOfInterest,     // Export only the area inside your polygon
  scale: 30,                  // Landsat resolution is 30 meters
  crs: 'EPSG:4326',           // Standard Coordinate Reference System (WGS84)
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true      // Makes the file easier to open in QGIS/ArcGIS
  }
});

 

Landsat 7

// 1. Define the Study Area
var areaOfInterest = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[10.904865172222156, 43.01615135702198],
          [10.904865172222156, 43.00648615894694],
          [10.925807860210437, 43.00648615894694],
          [10.925807860210437, 43.01615135702198]]], null, false);
// 2. Pre-processing Function for Landsat 8
var processL8 = function(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var ndvi = opticalBands.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');
  return image.addBands(ndvi).select('NDVI');
};

// 3. Load 2025 Median NDVI
var l7_2005 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
  .filterBounds(areaOfInterest)
  .filterDate('2005-01-01', '2005-12-31')
  .map(function(image) {
    var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
    // Landsat 7 NDVI: (B4 - B3) / (B4 + B3)
    var ndvi = opticalBands.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI');
    return image.addBands(ndvi).select('NDVI');
  })
  .median() // Create a clean composite for the year
  .clip(areaOfInterest);

// 4. Visualization Parameters
var palette = [
  '#FFFFFF', '#CE7E45', '#FCD163', '#66A000', '#207401', '#056201', '#011301'
];
var ndviViz = {min: 0, max: 0.9, palette: palette};

// 5. Display the Map
Map.addLayer(l7_2005, ndviViz, 'NDVI Median 2005');

// --- 6. ADDING THE LEGEND ---

// Create the legend panel
var legend = ui.Panel({
  style: {
    position: 'bottom-left',
    padding: '8px 15px'
  }
});

// Create legend title
var legendTitle = ui.Label({
  value: 'NDVI Index (2025)',
  style: {fontWeight: 'bold', fontSize: '16px', margin: '0 0 4px 0', padding: '0'}
});
legend.add(legendTitle);

// Helper function to create rows in the legend
var makeRow = function(color, name) {
  var colorBox = ui.Label({
    style: {
      backgroundColor: color,
      padding: '8px',
      margin: '0 0 4px 0'
    }
  });
  var description = ui.Label({
    value: name,
    style: {margin: '0 0 4px 6px'}
  });
  return ui.Panel({
    widgets: [colorBox, description],
    layout: ui.Panel.Layout.Flow('horizontal')
  });
};

// Define labels for the palette colors
var names = ['Low/No Veg (< 0.1)', 'Soil / Sparse', 'Moderate Veg', 'Dense Veg', 'Healthy Forest', 'Peak Vigour (> 0.8)'];
var colors = ['#FFFFFF', '#CE7E45', '#FCD163', '#66A000', '#207401', '#011301'];

// Add rows to legend
for (var i = 0; i < 6; i++) {
  legend.add(makeRow(colors[i], names[i]));
}

// Add legend to the map
Map.add(legend);

Export.image.toDrive({
  image: l7_2005,                // Updated to the 2005 variable
  description: 'NDVI_L7_2005_Tuscany', 
  folder: 'EarthEngine_Exports', 
  fileNamePrefix: 'NDVI_Landsat7_2005',
  region: areaOfInterest,     
  scale: 30,                  
  crs: 'EPSG:4326',           
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true      
  }
});


 

 

Analisi multitemporale NDVI con Google Earth Engine

Qualche esperimento con Google Earth Engine per vedere come si modifica nel tempo NDVI (arco di decenni fornito da Landsat)

 


 

// 1. Coordinate e Geometria
//var lon = 10.913507726271437;
//var lat = 43.01194709215584;

var lon = 10.926468160041452;
var lat = 43.01919548868449;


var poi = ee.Geometry.Point([lon, lat]);

// 2. Funzione di pre-elaborazione (Mask Cloud e Scaling)
// Valida per Landsat 5 e Landsat 7 (Collection 2 Level 2)
function preprocessLandsat(image) {
  var qa = image.select('QA_PIXEL');
  
  // Maschera per nuvole, ombre e cirri
  var mask = qa.bitwiseAnd(1 << 1).eq(0)
    .and(qa.bitwiseAnd(1 << 2).eq(0))
    .and(qa.bitwiseAnd(1 << 3).eq(0))
    .and(qa.bitwiseAnd(1 << 4).eq(0));

  // Fattori di scala per Collection 2
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  
  // Calcolo NDVI (NIR è B4, RED è B3 per entrambi i sensori L5 e L7)
  var ndvi = opticalBands.normalizedDifference(['SR_B4', 'SR_B3']).rename('NDVI');
  
  return image.addBands(ndvi).updateMask(mask);
}

// 3. Caricamento Collezione LANDSAT 5 (1985 - 2010)
var l5Col = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2")
  .filterBounds(poi)
  .filterDate('1985-01-01', '2010-01-01')
  .map(preprocessLandsat);

// 4. Caricamento Collezione LANDSAT 7 (2010 - 2024)
var l7Col = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2")
  .filterBounds(poi)
  .filterDate('2010-01-01', '2024-01-01')
  .map(preprocessLandsat);

// 5. Unione delle due collezioni
var mergedCol = l5Col.merge(l7Col);

// 6. Creazione del Grafico Temporale
var chart = ui.Chart.image.series({
  imageCollection: mergedCol.select('NDVI'),
  region: poi,
  reducer: ee.Reducer.mean(),
  scale: 30,
  xProperty: 'system:time_start'
})
.setOptions({
  title: 'Serie Storica NDVI: Landsat 5 e Landsat 7 (1985-2024)',
  vAxis: {title: 'NDVI', viewWindow: {min: -0.1, max: 1}},
  hAxis: {title: 'Anno', format: 'YYYY'},
  lineWidth: 1,
  pointSize: 2,
  trendlines: {0: {color: 'red', lineWidth: 1, opacity: 0.7}}, // Linea di tendenza globale
  series: {
    0: {color: '228B22'} // Colore verde foresta
  }
});

// 7. Visualizzazione
print(chart);
Map.centerObject(poi, 13);
Map.addLayer(poi, {color: 'red'}, 'Sito di analisi'); 

 

 

uno script piu' evoluto si puo'trovare qui. L'interesse e' che non viene calcolato un singolo punto ma il valore medio di un'area che risulta sicuramente piu' significativo
 

modificando per l'area di interesse

 

// 1. Definizione dell'Area di Studio (Rettangolo)
var areaOfInterest = ee.Geometry.Polygon([
  [[10.910535746643099, 43.01138938313477],
   [10.91448395831302, 43.01138938313477],
   [10.91448395831302, 43.00825125524326],
   [10.910535746643099, 43.00825125524326],
   [10.910535746643099, 43.01138938313477]]
]);

Map.centerObject(areaOfInterest, 15);
Map.addLayer(areaOfInterest, {color: 'd63000'}, 'Area di Studio');

// 2. Funzione Processamento con Soglia NDVI > 0.55
var processLandsat = function(image, nir, red) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var ndvi = opticalBands.normalizedDifference([nir, red]).rename('NDVI');
  
  // FILTRO: Manteniamo solo i pixel dove NDVI è maggiore di 0.55
  var maskThreshold = ndvi.gt(0.55);
  
  return image.addBands(ndvi).select('NDVI')
    .updateMask(maskThreshold) // Rimuove i valori bassi
    .copyProperties(image, ['system:time_start']);
};

// 3. Caricamento Collezioni
var L5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2").filterBounds(areaOfInterest)
  .map(function(img) { return processLandsat(img, 'SR_B4', 'SR_B3'); });

var L7 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2").filterBounds(areaOfInterest)
  .map(function(img) { return processLandsat(img, 'SR_B4', 'SR_B3'); });

var L8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").filterBounds(areaOfInterest)
  .map(function(img) { return processLandsat(img, 'SR_B5', 'SR_B4'); });

var L9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2").filterBounds(areaOfInterest)
  .map(function(img) { return processLandsat(img, 'SR_B5', 'SR_B4'); });

var mergedNDVI = L5.merge(L7).merge(L8).merge(L9);

// 4. Composizione Mensile (Massimo Valore)
var years = ee.List.sequence(1984, 2024);
var months = ee.List.sequence(1, 12);

var monthlyNDVI = ee.ImageCollection.fromImages(
  years.map(function(y) {
    return months.map(function(m) {
      var d = ee.Date.fromYMD(y, m, 1);
      return mergedNDVI
        .filter(ee.Filter.calendarRange(y, y, 'year'))
        .filter(ee.Filter.calendarRange(m, m, 'month'))
        .max() // Seleziona il valore più alto (più verde) del mese
        .set('system:time_start', d.millis());
    });
  }).flatten()
);

// Rimuoviamo i periodi senza dati (es. inverno o periodi sotto soglia)
var monthlyNDVI_Clean = monthlyNDVI.filter(ee.Filter.listContains('system:band_names', 'NDVI'));

// 5. Grafico
var chart = ui.Chart.image.series({
  imageCollection: monthlyNDVI_Clean,
  region: areaOfInterest,
  reducer: ee.Reducer.mean(),
  scale: 30,
  xProperty: 'system:time_start'
})
.setOptions({
  title: 'NDVI Filtrato (>0.55) - Peak Vegetation Trend',
  vAxis: {title: 'NDVI', viewWindow: {min: 0.55, max: 0.9}},
  hAxis: {title: 'Anno', format: 'YYYY'},
  trendlines: {0: {color: 'red', lineWidth: 1}},
  series: {0: {color: '228B22', pointSize: 2, lineWidth: 1}}
});

print(chart); 

 


 


in generale sembra che NDVI sia migliorato a partire dal 2013

domenica 8 febbraio 2026

Flauto traverso barocco stampato 3D

Ho trovato questo modello di fluato traverso barocco da stampare https://www.thingiverse.com/thing:4572430 che e' una copia di un modello del 1740 di Bizet

  


I tenoni sono abbastanza sprecisi...ero partito per usare il nastro di teflon per garantire la tenuta ma alla fine per fare spessore ho usato in nastro carta. Per il tampone della chiave del Reb ho usato nastro biadesivo 

 

 



Io non saro' certo un grande flautista ma io non riesco a suonarlo. Non c'e' volume 

Inoltre la distanza dei fori e' diabolicamente molto superiore a quella di un flauto moderno. Peraltro la diteggiatura barocca, spesso con forchette, non aiuta 


 Esiste un altro modello disponibile da stampare che e' una copia di un flauto del 1750 conservato ad Amsterdam ed opera di Deppe  https://www.rijksmuseum.nl/en/collection/object/Flute--6e8161fcef916a7207f25a0812507870

Di questo secondo flauto esiste un video con nei commenti il link per i file STL.  

 


 

 

 

lunedì 2 febbraio 2026

Correzione atmosferica su dati EnMap

 Mi e' stato girato un articolo  

Detecting rare earth elements using EnMAP hyperspectral satellite data: a case study from Mountain Pass, California

in cui si indica che e' possibile mappare il neodimio da satellite usando bande VNIR. Mi sono incuriosito ed ho provato a scaricare i dati pubblici di Enmap relativi alla piu' grande miniera di Bayan Obo in Cina

Longitude:    109°54'55" E    degree
Latitude:    41°40'59" N    degree

 Dal punto di vista degli assorbimenti si vede che gli assorbimenti significativi sono a circa 740, 798 ed 870 nm (per rimanere nel campo del VNIR)

Ho provato a selezionare alcuni punti sia all'interno dell'area di escavazione che all'esterno
 
 

e questi sono gli spettri


 come si vede c'e' un picco di assorbimento a 760 nm ma questo e' presente sia nei dati dentro alla miniera che all'esterno ed i picchi sono via via che lo spettro ha mediamente una riflettanza maggiore

piu' in dettaglio 



 come si vede nella zona di 760 c'e' un profondo picco di assorbimento dell'ossigeno.

I dati che ho scaricato da Enmap erano di livello L2A quindi gia' trattati per la rimozione del contributo atmosferico  

L'ipotesi che l'algoritmo di rimozione atmosferica non abbia eliminato completamente il contributo dell'ossigeno e che quindi i dati non siano significativi nell'intorno della lunghezza d'onda di 760 nm 

 

 

 


domenica 1 febbraio 2026

Covox Speech Thing

Nel cassonetto delle pile usate a lavoro ho trovato questo circuito ed e' stato un tutto negli anni 80

Non lo ho mai usato (sono partito da un clone Soundblaster ad 8 bit) ma questa e' una primordiale scheda audio da collegare alla porta parallela

  


Originariamente nota come Covox Speech Thing (questa e' un clone) e' di fatto un convertitore Digitale Analogico ad 8 Bit basato solo resistenze in configurazione R-2R (https://en.wikipedia.org/wiki/Resistor_ladder)

Al termine c'e' uno stadio di amplificazione (da qui la necessita' delle pile) e l'uscita audio..nessun integrato..talmente semplice ed economico che se non ricordo male era incluso nelle confezioni di alcuni giochi

 Il problema di questi aggeggi e' che non avendo un timer se il processore del computer era impiegato in altri compiti oltre a quelli di aggiornare il suono la musica rallentava 

 

 

 



 

 

 

 

 

 

venerdì 23 gennaio 2026

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()

mercoledì 14 gennaio 2026

Microfoni Kinect1

Ho trovato quasi per caso su Internet che il Kinect 1 ha un array di 4 microfoni che permettono di determinare la direzione del suono (ed essendo tutto integrato non ha problemi di jitter) 

Bus 001 Device 016: ID 045e:02b0 Microsoft Corp. Xbox NUI Motor
Bus 001 Device 018: ID 045e:02bb Microsoft Corp. Kinect Audio
Bus 001 Device 019: ID 045e:02ae Microsoft Corp. Xbox NUI Camera
Bus 002 Device 001: ID 1d6b:0003 Linux Foundation 3.0 root hub

L'array di 4 microfoni e' indicato da 0 (lato destro) a 3 (lato sinistro) con una distanza tra 0 e 3 di 22.6 cm. I microfoni 0,1,2 sono molto ravvicinati tra di loro, solo il 3 e' distanziato alla sinistra


 

Una delle limitazione di kinect e' che il campionamento massimo e' di 16k il che limita la precisione a 2.14 cm ed 24 bit (anche se sono trasmessi come 32 bit)

arecord -l

 card 1: Audio [Kinect USB Audio], device 0: USB Audio [USB Audio]
  Subdevices: 1/1
  Subdevice #0: subdevice #0


arecord -D hw:1,0 -c 4 -r 16000 -f S32_LE -t raw kinect_array.raw 

ffmpeg -f s32le -ar 16000 -ac 4 -i kinect_array.raw kinect_output.wav 

 

import numpy as np
import subprocess
import math
from scipy.signal import correlate

# --- CONFIGURAZIONE ---
CHANNELS = 4
RATE = 16000
CHUNK_SIZE = 1024 * 4 # Dimensione del blocco da analizzare (circa 0.25 secondi)
v = 343 # Velocità del suono m/s
d = 0.226 # Distanza Mic 0 - Mic 3 (metri)

# Comando per avviare arecord e inviare i dati RAW alla stdout
cmd = [
'arecord', '-D', 'hw:1,0', '-c', str(CHANNELS), '-r', str(RATE),
'-f', 'S32_LE', '-t', 'raw', '-q'
]

# Avvio del processo
proc = subprocess.Popen(cmd, stdout=subprocess.PIPE)

print("Ascolto in corso... Muoviti davanti al Kinect o schiocca le dita!")
print("Premi Ctrl+C per fermare.\n")

try:
while True:
# Leggi un blocco di dati (4 canali * 4 byte per campione * CHUNK_SIZE)
raw_bytes = proc.stdout.read(CHUNK_SIZE * CHANNELS * 4)
if not raw_bytes:
break

# Converti byte in array numpy
data = np.frombuffer(raw_bytes, dtype=np.int32).reshape(-1, CHANNELS)

# Prendi Mic 0 e Mic 3
mic_0 = data[:, 0].astype(np.float32)
mic_3 = data[:, 3].astype(np.float32)

# Calcola correlazione
corr = correlate(mic_0, mic_3, mode='full')
lag = np.argmax(corr) - (len(mic_0) - 1)
# Calcolo tempo e angolo
dt = lag / RATE
# Limita l'argomento di arcsin tra -1 e 1 per evitare errori
sin_argument = (v * dt) / d
if abs(sin_argument) <= 1.0:
angle_rad = math.asin(sin_argument)
angle_deg = math.degrees(angle_rad)
# Visualizzazione a barre dinamica
bar_position = int(angle_deg / 10) + 9
display = ["-"] * 19
display[bar_position] = "O"
print(f"[{''.join(display)}] Angolo: {angle_deg:6.1f}°", end="\r")

except KeyboardInterrupt:
print("\nFermato dall'utente.")
finally:
proc.terminate()

 

Mappa NDVI medio annuale con Google Earth Engine

Al contrario del precedente post https://debiaonoldcomputers.blogspot.com/2026/02/analisi-multitemporale-ndvi-con-google.html  qui l'ide...