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

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