Edge AI 기반 실시간 진동 분석: Raspberry Pi와 MEMS 센서로 구현하는 저비용 예지보전 시스템

# 필수 라이브러리
import adafruit_adxl34x
import board
import numpy as np
from scipy import signal
import tflite_runtime.interpreter as tflite

import busio
import adafruit_adxl34x

# I2C 버스 초기화
i2c = busio.I2C(board.SCL, board.SDA)
accel = adafruit_adxl34x.ADXL345(i2c)

# 샘플링 레이트 설정 (1600Hz)
accel.data_rate = adafruit_adxl34x.DataRate.RATE_1600_HZ
accel.range = adafruit_adxl34x.Range.RANGE_16_G

import time

SAMPLE_RATE = 1600  # Hz
WINDOW_SIZE = 1024  # 샘플 수 (약 0.64초)

def collect_vibration_data():
    """1024개 샘플 수집 (약 0.64초 윈도우)"""
    buffer = []
    for _ in range(WINDOW_SIZE):
        x, y, z = accel.acceleration
        # 3축 벡터 크기 계산
        magnitude = np.sqrt(x**2 + y**2 + z**2)
        buffer.append(magnitude)
        time.sleep(1 / SAMPLE_RATE)
    return np.array(buffer)

def extract_frequency_features(signal_data, sample_rate=1600):
    """FFT를 통한 주파수 특징 추출"""
    # Hanning 윈도우 적용 (스펙트럼 누수 방지)
    windowed = signal_data * np.hanning(len(signal_data))

    # FFT 계산
    fft_vals = np.fft.rfft(windowed)
    fft_freq = np.fft.rfftfreq(len(windowed), 1/sample_rate)
    fft_mag = np.abs(fft_vals)

    # 주요 특징 추출
    features = {
        'peak_freq': fft_freq[np.argmax(fft_mag)],  # 최대 진폭 주파수
        'rms': np.sqrt(np.mean(signal_data**2)),     # RMS (진동 강도)
        'kurtosis': np.mean((signal_data - np.mean(signal_data))**4) / np.std(signal_data)**4,
        'spectral_centroid': np.sum(fft_freq * fft_mag) / np.sum(fft_mag)
    }
    return features, fft_mag[:128]  # 128개 주파수 빈 사용

import tensorflow as tf

def build_model():
    model = tf.keras.Sequential([
        tf.keras.layers.Input(shape=(128, 1)),  # 128개 주파수 빈
        tf.keras.layers.Conv1D(16, 5, activation='relu'),
        tf.keras.layers.MaxPooling1D(2),
        tf.keras.layers.Conv1D(32, 5, activation='relu'),
        tf.keras.layers.GlobalAveragePooling1D(),
        tf.keras.layers.Dense(16, activation='relu'),
        tf.keras.layers.Dense(1, activation='sigmoid')  # 이상/정상 (0/1)
    ])
    return model

model = build_model()
model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])

# 양자화(Quantization)로 모델 크기 1/4 축소
converter = tf.lite.TFLiteConverter.from_keras_model(model)
converter.optimizations = [tf.lite.Optimize.DEFAULT]
tflite_model = converter.convert()

with open('vibration_model.tflite', 'wb') as f:
    f.write(tflite_model)

import tflite_runtime.interpreter as tflite

# 모델 로드
interpreter = tflite.Interpreter(model_path='vibration_model.tflite')
interpreter.allocate_tensors()

input_details = interpreter.get_input_details()
output_details = interpreter.get_output_details()

def predict_anomaly(fft_spectrum):
    """실시간 이상 탐지"""
    # 입력 데이터 정규화
    input_data = fft_spectrum.reshape(1, 128, 1).astype(np.float32)
    input_data = (input_data - np.mean(input_data)) / np.std(input_data)

    # 추론 실행
    interpreter.set_tensor(input_details[0]['index'], input_data)
    interpreter.invoke()
    output = interpreter.get_tensor(output_details[0]['index'])

    return output[0][0]  # 이상 확률 (0~1)

def calculate_bpfo(rpm, n_balls, ball_dia, pitch_dia, contact_angle=0):
    """외륜 결함 주파수 계산"""
    fs = rpm / 60  # Hz 변환
    return (n_balls / 2) * fs * (1 - (ball_dia / pitch_dia) * np.cos(contact_angle))

# 예시: 1800 RPM 모터, 9개 볼 베어링
bpfo = calculate_bpfo(1800, 9, 7.9, 39.0)
print(f"외륜 결함 주파수: {bpfo:.2f} Hz")  # 약 120 Hz

import paho.mqtt.client as mqtt

def send_alert(anomaly_score, features):
    """MQTT로 경보 전송"""
    client = mqtt.Client()
    client.connect("mqtt.local", 1883)

    payload = {
        "timestamp": time.time(),
        "anomaly_score": float(anomaly_score),
        "peak_freq": features['peak_freq'],
        "rms": features['rms'],
        "kurtosis": features['kurtosis']
    }

    client.publish("factory/vibration/alert", json.dumps(payload))
    client.disconnect()

# 메인 루프
while True:
    signal_data = collect_vibration_data()
    features, fft_spectrum = extract_frequency_features(signal_data)
    anomaly_score = predict_anomaly(fft_spectrum)

    if anomaly_score > 0.75:  # 임계값 75%
        send_alert(anomaly_score, features)
        print(f"⚠️  이상 탐지! (확률: {anomaly_score:.2%})")

    time.sleep(10)  # 10초마다 분석

# 적응형 샘플링: 정상 시 10초, 이상 시 1초 간격
if anomaly_score < 0.3:
    sleep_time = 10
elif anomaly_score < 0.7:
    sleep_time = 5
else:
    sleep_time = 1  # 고위험 시 실시간 모니터링

# 온도 센서 추가 (DS18B20)
import w1thermsensor

temp_sensor = w1thermsensor.W1ThermSensor()
temperature = temp_sensor.get_temperature()

# 진동 + 온도 융합 분석
if anomaly_score > 0.7 and temperature > 70:
    alert_level = "CRITICAL"  # 베어링 과열 + 진동 이상