English
sofia-engine
edge-ai
industrial-ai
scientific-computing
embedded-ai
signal-processing
digital-signal-processing
predictive-maintenance
condition-monitoring
vibration-analysis
anomaly-detection
industrial-iot
iiot
telemetry
edge-computing
tinyml
on-device-learning
embedded-systems
machine-health
time-series
python
typescript
c
File size: 4,063 Bytes
876458a | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 | """Scientific Signal Analysis Example for Sofia Engine.
Demonstrates deterministic digital signal processing across mechanical vibration
and industrial electrical domains:
1. FFT magnitude spectrum & peak detection
2. Welch power spectral density & Parseval energy conservation
3. Hilbert transform analytic amplitude envelope
4. Fortescue 3-phase symmetrical components & Voltage Unbalance Factor (VUF)
Run:
python examples/signal_analysis.py
"""
from __future__ import annotations
import numpy as np
from sofia_ai.signal.electrical import compute_symmetrical_components
from sofia_ai.signal.envelope import amplitude_envelope
from sofia_ai.signal.spectral import peak_frequency, rfft_magnitude, welch_psd
def run_signal_analysis() -> None:
print("=== Sofia Engine: Scientific Signal Analysis ===")
# -------------------------------------------------------------------------
# 1. Vibration Spectral Analysis (FFT & Welch PSD)
# -------------------------------------------------------------------------
fs = 2000.0 # 2 kHz sampling rate
n = 2000 # 1 second duration
t = np.arange(n) / fs
# Synthesize multi-tone machine vibration:
# 50 Hz shaft rotation (0.8 m/s^2) + 250 Hz gear-mesh tone (0.3 m/s^2)
tone_50 = 0.8 * np.sin(2 * np.pi * 50.0 * t)
tone_250 = 0.3 * np.sin(2 * np.pi * 250.0 * t)
noise = 0.05 * np.random.default_rng(123).normal(size=n)
signal = tone_50 + tone_250 + noise
print(f"\n[1] Vibration Telemetry ({n} samples @ {fs} Hz)")
# Compute FFT magnitude spectrum
spec = rfft_magnitude(signal, fs, window="hann", detrend_mean=True)
f_peak = peak_frequency(spec.frequencies, spec.values)
mag_peak = float(np.max(spec.values))
print(f" FFT Peak Frequency: {f_peak:.1f} Hz (Magnitude: {mag_peak:.4f})")
# Compute Welch PSD
psd = welch_psd(signal, fs, segment_length=512, overlap=0.5, window="hann")
# Parseval check: integral of PSD approx equals variance of signal
signal_variance = float(np.var(signal))
psd_total_power = psd.total_power
parseval_ratio = psd_total_power / max(signal_variance, 1e-12)
print(f" Welch PSD Total Power: {psd_total_power:.4f}")
print(f" Time-Domain Variance: {signal_variance:.4f}")
print(f" Parseval Energy Ratio: {parseval_ratio:.4f} (expected ~1.0)")
# -------------------------------------------------------------------------
# 2. Bearing Fault Envelope Demodulation (Hilbert Transform)
# -------------------------------------------------------------------------
print("\n[2] Amplitude Demodulation (Hilbert Analytic Envelope)")
# Amplitude modulated signal: 200 Hz carrier modulated at 15 Hz fault pass frequency
carrier_freq = 200.0
mod_freq = 15.0
envelope_true = 1.0 + 0.6 * np.cos(2 * np.pi * mod_freq * t)
am_signal = envelope_true * np.cos(2 * np.pi * carrier_freq * t)
env = amplitude_envelope(am_signal)
env_spec = rfft_magnitude(env - np.mean(env), fs, window="hann")
env_peak_f = peak_frequency(env_spec.frequencies, env_spec.values)
print(f" Extracted Modulation Frequency: {env_peak_f:.1f} Hz (expected {mod_freq:.1f} Hz)")
# -------------------------------------------------------------------------
# 3. 3-Phase Symmetrical Components (Fortescue Transformation)
# -------------------------------------------------------------------------
print("\n[3] 3-Phase Symmetrical Components (IEC / IEEE 519 analysis)")
# Unbalanced 3-phase voltages (Va=230V @ 0°, Vb=220V @ -125°, Vc=215V @ 115°)
sym = compute_symmetrical_components(
va_amp=230.0, va_phase_deg=0.0,
vb_amp=220.0, vb_phase_deg=-125.0,
vc_amp=215.0, vc_phase_deg=115.0,
)
print(f" Positive Sequence (V1): {sym.v1_pos_seq_v:.2f} V")
print(f" Negative Sequence (V2): {sym.v2_neg_seq_v:.2f} V")
print(f" Zero Sequence (V0): {sym.v0_zero_seq_v:.2f} V")
print(f" Voltage Unbalance Factor (VUF): {sym.vuf_percent:.2f}%")
if __name__ == "__main__":
run_signal_analysis()
|