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