SciPy: Scientific Computing dengan Python

1. Pengenalan SciPy

SciPy (Scientific Python) adalah library Python yang menyediakan algoritma numerik dan fungsi matematika untuk scientific computing. Dibangun di atas NumPy, SciPy melengkapi ekosistem scientific Python dengan modul-modul khusus untuk optimization, integration, interpolation, signal processing, linear algebra, dan banyak lagi.

SciPy adalah fondasi dari banyak riset sains dan engineering. Mulai dari analisis data eksperimen fisika, simulasi struktur jembatan, processing sinyal audio, hingga optimasi portofolio keuangan — semua bisa dilakukan dengan SciPy.

Modul SciPy

Modul	Fungsi	Contoh Penggunaan
scipy.optimize	Optimasi & root finding	Minimalkan biaya, cari akar persamaan
scipy.interpolate	Interpolasi data	Isi data yang hilang, kurva fitting
scipy.signal	Signal processing	Filter audio, FFT, deteksi puncak
scipy.linalg	Linear algebra	Matrix decomposition, eigenvalues
scipy.stats	Statistik	Distribusi probabilitas, uji hipotesis
scipy.integrate	Numerical integration	Integral, ODE solver
scipy.sparse	Sparse matrices	Graph problems, large-scale computation
scipy.spatial	Spatial algorithms	Nearest neighbor, convex hull, KDTree
scipy.ndimage	N-dimensional image	Image processing, filtering

# Instalasi
pip install scipy numpy matplotlib

# Import modul yang dibutuhkan
import numpy as np
from scipy import optimize, interpolate, signal, linalg, stats, integrate
import matplotlib.pyplot as plt

2. Optimization

Modul scipy.optimize menyediakan berbagai algoritma untuk menemukan minimum/maximum fungsi, mencari akar persamaan, dan curve fitting. Ini sangat penting dalam machine learning, engineering, dan riset operasi.

Minimisasi Fungsi

import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# Fungsi yang akan diminimalkan: f(x) = (x-3)^2 + 2*x + 5
def objective(x):
    return (x - 3)**2 + 2*x + 5

# Minimisasi dari tebakan awal x=0
hasil = minimize(objective, x0=0, method='BFGS')
print(f"Minimum di x = {hasil.x[0]:.4f}")
print(f"Nilai minimum = {hasil.fun:.4f}")
print(f"Konvergen? {hasil.success}")
print(f"Iterasi: {hasil.nit}")

# Visualisasi
x = np.linspace(-2, 8, 100)
y = [objective(xi) for xi in x]
plt.figure(figsize=(10, 6))
plt.plot(x, y, 'b-', linewidth=2, label='f(x) = (x-3)² + 2x + 5')
plt.plot(hasil.x[0], hasil.fun, 'r*', markersize=15, label=f'Minimum ({hasil.x[0]:.2f}, {hasil.fun:.2f})')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Optimasi Fungsi dengan SciPy')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('optimization.png', dpi=150)
plt.show()

Optimasi Multivariabel

import numpy as np
from scipy.optimize import minimize

# Fungsi Rosenbrock: benchmark standar untuk optimasi
# f(x,y) = (1-x)^2 + 100*(y-x^2)^2
def rosenbrock(params):
    x, y = params
    return (1 - x)**2 + 100 * (y - x**2)**2

# Minimisasi dari titik awal
hasil = minimize(rosenbrock, x0=[0, 0], method='Nelder-Mead')
print(f"Minimum di x={hasil.x[0]:.6f}, y={hasil.x[1]:.6f}")
print(f"Nilai minimum = {hasil.fun:.10f}")  # Harusnya ~0 di (1,1)

# Metode L-BFGS-B dengan batasan (bounds)
from scipy.optimize import minimize, Bounds

bounds = Bounds([0, 0], [5, 5])  # x,y antara 0 dan 5
hasil = minimize(rosenbrock, x0=[0, 0], method='L-BFGS-B', bounds=bounds)
print(f"Dengan batasan: x={hasil.x[0]:.4f}, y={hasil.x[1]:.4f}")

Root Finding

import numpy as np
from scipy.optimize import root, root_scalar

# Cari akar: x^3 - 6x^2 + 11x - 6 = 0
# (Akar: x=1, x=2, x=3)

def equation(x):
    return x**3 - 6*x**2 + 11*x - 6

# Single variable root
for x0 in [0, 1.5, 2.5, 4]:
    sol = root_scalar(equation, bracket=[x0-0.5, x0+0.5], method='brentq')
    print(f"  Akar near {x0}: x = {sol.root:.6f}")

# Multivariable root finding
def equations(vars):
    x, y = vars
    eq1 = x**2 + y**2 - 4  # x^2 + y^2 = 4
    eq2 = x - y              # x = y
    return [eq1, eq2]

sol = root(equations, x0=[1, 1])
print(f"\nSolusi: x={sol.x[0]:.4f}, y={sol.x[1]:.4f}")
print(f"Verifikasi: {sol.fun}")  # Harusnya ~[0, 0]

Curve Fitting

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Data eksperimen (suhu vs waktu)
waktu = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
suhu = np.array([95, 82, 71, 62, 55, 49, 44, 40, 37, 35, 34])

# Model: pendinginan Newton T(t) = T_env + (T0 - T_env) * exp(-k*t)
def model(t, T_env, T0, k):
    return T_env + (T0 - T_env) * np.exp(-k * t)

# Fit parameter
popt, pcov = curve_fit(model, waktu, suhu, p0=[30, 95, 0.3])
T_env, T0, k = popt
perr = np.sqrt(np.diag(pcov))  # Standard error

print(f"T_env = {T_env:.2f} ± {perr[0]:.2f}")
print(f"T0 = {T0:.2f} ± {perr[1]:.2f}")
print(f"k = {k:.4f} ± {perr[2]:.4f}")

# Visualisasi
t_fit = np.linspace(0, 15, 100)
plt.figure(figsize=(10, 6))
plt.scatter(waktu, suhu, color='red', label='Data Eksperimen')
plt.plot(t_fit, model(t_fit, *popt), 'b-', label=f'Fit: T={T_env:.1f}+{T0-T_env:.1f}*e^(-{k:.3f}t)')
plt.xlabel('Waktu (menit)')
plt.ylabel('Suhu (°C)')
plt.title('Pendinginan Newton - Curve Fitting')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('curve_fitting.png', dpi=150)
plt.show()

3. Interpolation

Interpolasi digunakan untuk memperkirakan nilai di antara titik-titik data yang diketahui. Sangat berguna untuk mengisi data yang hilang, resampling, atau membuat kurva halus dari data kasar.

import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt

# Data asli (kasar)
x_data = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8])
y_data = np.array([0, 0.8, 0.9, 0.1, -0.8, -1.0, -0.3, 0.7, 1.0])

# Berbagai metode interpolasi
x_fine = np.linspace(0, 8, 200)

# 1. Linear Interpolation
f_linear = interpolate.interp1d(x_data, y_data, kind='linear')
y_linear = f_linear(x_fine)

# 2. Cubic Spline
f_cubic = interpolate.interp1d(x_data, y_data, kind='cubic')
y_cubic = f_cubic(x_fine)

# 3. Akima (menghindari overshooting)
f_akima = interpolate.Akima1DInterpolator(x_data, y_data)
y_akima = f_akima(x_fine)

# Visualisasi
plt.figure(figsize=(12, 6))
plt.scatter(x_data, y_data, color='red', s=100, zorder=5, label='Data Asli')
plt.plot(x_fine, y_linear, '--', label='Linear', alpha=0.7)
plt.plot(x_fine, y_cubic, '-', linewidth=2, label='Cubic Spline')
plt.plot(x_fine, y_akima, '-.', label='Akima', alpha=0.7)
plt.legend()
plt.title('Perbandingan Metode Interpolasi')
plt.grid(True, alpha=0.3)
plt.savefig('interpolation.png', dpi=150)
plt.show()

# Interpolasi 2D
from scipy.interpolate import griddata

# Data irregular 2D
points = np.random.rand(100, 2)
values = np.sin(points[:, 0] * 2 * np.pi) * np.cos(points[:, 1] * 2 * np.pi)

# Grid reguler
grid_x, grid_y = np.mgrid[0:1:100j, 0:1:100j]
grid_z = griddata(points, values, (grid_x, grid_y), method='cubic')

print(f"Interpolasi 2D: grid shape = {grid_z.shape}")

4. Signal Processing

Modul scipy.signal adalah toolkit lengkap untuk analisis dan pemrosesan sinyal — dari filter digital hingga transformasi Fourier.

FFT (Fast Fourier Transform)

import numpy as np
from scipy import signal, fft
import matplotlib.pyplot as plt

# Buat sinyal: gabungan 3 frekuensi + noise
fs = 1000  # Sample rate (Hz)
t = np.arange(0, 1.0, 1/fs)

# Sinyal: 50Hz + 120Hz + 200Hz + noise
sinyal = (0.7 * np.sin(2*np.pi*50*t) +
          1.0 * np.sin(2*np.pi*120*t) +
          0.5 * np.sin(2*np.pi*200*t) +
          0.3 * np.random.randn(len(t)))

# FFT
N = len(sinyal)
yf = fft.fft(sinyal)
xf = fft.fftfreq(N, 1/fs)[:N//2]
amplitude = 2.0/N * np.abs(yf[:N//2])

# Visualisasi
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
ax1.plot(t[:200], sinyal[:200])
ax1.set_title('Domain Waktu')
ax1.set_xlabel('Waktu (s)')
ax1.set_ylabel('Amplitudo')

ax2.plot(xf, amplitude)
ax2.set_title('Domain Frekuensi (FFT)')
ax2.set_xlabel('Frekuensi (Hz)')
ax2.set_ylabel('Amplitudo')
ax2.set_xlim(0, 300)
plt.tight_layout()
plt.savefig('fft.png', dpi=150)
plt.show()

# Identifikasi frekuensi dominan
top_indices = np.argsort(amplitude)[-3:][::-1]
for idx in top_indices:
    print(f"  Frekuensi: {xf[idx]:.1f} Hz, Amplitudo: {amplitude[idx]:.3f}")

Filter Digital

import numpy as np
from scipy.signal import butter, filtfilt, iirnotch
import matplotlib.pyplot as plt

# Sinyal dengan noise frekuensi tinggi
fs = 1000
t = np.arange(0, 1.0, 1/fs)
sinyal_bersih = np.sin(2*np.pi*5*t)  # 5 Hz signal
noise = 0.5 * np.sin(2*np.pi*100*t) + 0.3 * np.random.randn(len(t))
sinyal_noisy = sinyal_bersih + noise

# Low-pass Butterworth filter (cutoff 20 Hz)
def lowpass_filter(data, cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return filtfilt(b, a, data)

sinyal_filtered = lowpass_filter(sinyal_noisy, cutoff=20, fs=fs)

# Notch filter (hilangkan 50Hz interference)
def notch_filter(data, freq, fs, Q=30):
    b, a = iirnotch(freq, Q, fs)
    return filtfilt(b, a, data)

sinyal_notch = notch_filter(sinyal_noisy, 100, fs)

# Visualisasi
fig, axes = plt.subplots(3, 1, figsize=(12, 10), sharex=True)
axes[0].plot(t, sinyal_noisy, 'gray', alpha=0.7)
axes[0].plot(t, sinyal_bersih, 'b-', linewidth=2)
axes[0].set_title('Sinyal Asli + Noise')

axes[1].plot(t, sinyal_filtered, 'r-', linewidth=2)
axes[1].set_title('Low-Pass Filter (cutoff 20Hz)')

axes[2].plot(t, sinyal_notch, 'g-', linewidth=2)
axes[2].set_title('Notch Filter (remove 100Hz)')

for ax in axes:
    ax.set_xlim(0, 0.5)
plt.tight_layout()
plt.savefig('filters.png', dpi=150)
plt.show()

Deteksi Puncak (Peak Detection)

import numpy as np
from scipy.signal import find_peaks, peak_prominences
import matplotlib.pyplot as plt

# Sinyal dengan beberapa puncak
t = np.linspace(0, 10, 1000)
sinyal = (np.sin(t) * np.exp(-0.1*t) +
          0.5*np.sin(3*t) * np.exp(-0.2*t) +
          0.3*np.random.randn(len(t)))

# Deteksi puncak
peaks, properties = find_peaks(sinyal,
    height=0.3,        # Min height
    distance=50,        # Min distance antar puncak
    prominence=0.2,     # Min prominence
)

# Hitung prominences
prominences = peak_prominences(sinyal, peaks)[0]

print(f"Ditemukan {len(peaks)} puncak:")
for i, (idx, prom) in enumerate(zip(peaks, prominences)):
    print(f"  Puncak {i+1}: t={t[idx]:.2f}s, amplitude={sinyal[idx]:.3f}, prominence={prom:.3f}")

plt.figure(figsize=(12, 5))
plt.plot(t, sinyal, 'b-', linewidth=1)
plt.plot(t[peaks], sinyal[peaks], 'rv', markersize=10, label='Puncak terdeteksi')
plt.xlabel('Waktu (s)')
plt.ylabel('Amplitudo')
plt.title('Peak Detection dengan SciPy')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('peaks.png', dpi=150)
plt.show()

5. Linear Algebra

scipy.linalg melengkapi NumPy dengan operasi linear algebra yang lebih lengkap dan seringkali lebih cepat.

import numpy as np
from scipy import linalg

# Matrix operations
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 0]])

# Determinan
det = linalg.det(A)
print(f"Determinan: {det:.0f}")  # 27

# Inverse
A_inv = linalg.inv(A)
print(f"A * A_inv = I? {np.allclose(A @ A_inv, np.eye(3))}")

# Eigenvalues & Eigenvectors
eigenvalues, eigenvectors = linalg.eig(A)
print(f"\nEigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")

# SVD (Singular Value Decomposition)
U, s, Vh = linalg.svd(A)
print(f"\nSingular values: {s}")

# Solve linear system: Ax = b
b = np.array([1, 2, 3])
x = linalg.solve(A, b)
print(f"\nSolusi Ax=b: x = {x}")
print(f"Verifikasi: A*x = {A @ x} (harus = {b})")

# LU Decomposition
P, L, U = linalg.lu(A)
print(f"\nLU Decomposition:")
print(f"P:\n{P}")
print(f"L:\n{L}")
print(f"U:\n{U}")
print(f"P@L@U = A? {np.allclose(P @ L @ U, A)}")

Matrix Decomposition Lanjutan

import numpy as np
from scipy import linalg

# Cholesky decomposition (untuk symmetric positive definite matrix)
A_sym = np.array([[4, 2, 1],
                  [2, 5, 3],
                  [1, 3, 6]])

L = linalg.cholesky(A_sym, lower=True)
print(f"Cholesky: A = L * L^T")
print(f"Verifikasi: {np.allclose(L @ L.T, A_sym)}")

# QR Decomposition
A = np.random.randn(5, 3)
Q, R = linalg.qr(A)
print(f"\nQR: Q shape={Q.shape}, R shape={R.shape}")
print(f"Q orthogonal? {np.allclose(Q.T @ Q, np.eye(3))}")

# Matrix norm
A = np.array([[1, -1], [2, 3]])
print(f"\nFrobenius norm: {linalg.norm(A, 'fro'):.4f}")
print(f"Spectral norm: {linalg.norm(A, 2):.4f}")

6. Statistik

scipy.stats menyediakan lebih dari 100 distribusi probabilitas dan berbagai uji statistik.

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

# Distribusi Normal
dist = stats.norm(loc=170, scale=10)  # Mean=170cm, StdDev=10cm

# Probabilitas
print(f"P(X < 160) = {dist.cdf(160):.4f}")  # CDF
print(f"P(X > 180) = {1 - dist.cdf(180):.4f}")
print(f"P(160 < X < 180) = {dist.cdf(180) - dist.cdf(160):.4f}")

# Percentiles
print(f"\nMedian: {dist.median():.1f} cm")
print(f"95th percentile: {dist.ppf(0.95):.1f} cm")

# Random sampling
samples = dist.rvs(size=1000, random_state=42)
print(f"Sample mean: {samples.mean():.2f}")
print(f"Sample std: {samples.std():.2f}")

# Visualisasi distribusi
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# Normal
x = np.linspace(130, 210, 200)
axes[0].plot(x, dist.pdf(x), 'b-', linewidth=2)
axes[0].fill_between(x, dist.pdf(x), where=(x < 160) | (x > 180), alpha=0.3)
axes[0].set_title('Distribusi Normal (μ=170, σ=10)')
axes[0].set_xlabel('Tinggi (cm)')

# Binomial
n, p = 20, 0.5
x_binom = np.arange(0, 21)
axes[1].bar(x_binom, stats.binom.pmf(x_binom, n, p), color='green', alpha=0.7)
axes[1].set_title(f'Distribusi Binomial (n={n}, p={p})')

# Poisson
lam = 5
x_poisson = np.arange(0, 15)
axes[2].bar(x_poisson, stats.poisson.pmf(x_poisson, lam), color='red', alpha=0.7)
axes[2].set_title(f'Distribusi Poisson (λ={lam})')

plt.tight_layout()
plt.savefig('distributions.png', dpi=150)
plt.show()

Hypothesis Testing

import numpy as np
from scipy import stats

np.random.seed(42)

# Data: kelompok kontrol vs treatment
kontrol = np.random.normal(loc=100, scale=15, size=50)
treatment = np.random.normal(loc=110, scale=15, size=50)

# 1. T-test (apakah rata-rata berbeda signifikan?)
t_stat, p_value = stats.ttest_ind(kontrol, treatment)
print(f"=== Independent T-Test ===")
print(f"t-statistic: {t_stat:.4f}")
print(f"p-value: {p_value:.6f}")
print(f"Signifikan (α=0.05)? {'Ya' if p_value < 0.05 else 'Tidak'}")

# 2. Uji normalitas (Shapiro-Wilk)
stat_k, p_k = stats.shapiro(kontrol)
stat_t, p_t = stats.shapiro(treatment)
print(f"\n=== Uji Normalitas (Shapiro-Wilk) ===")
print(f"Kontrol: W={stat_k:.4f}, p={p_k:.4f}")
print(f"Treatment: W={stat_t:.4f}, p={p_t:.4f}")

# 3. Chi-square test
observed = np.array([[50, 30, 20], [35, 40, 25]])
chi2, p_val, dof, expected = stats.chi2_contingency(observed)
print(f"\n=== Chi-Square Test ===")
print(f"Chi2 = {chi2:.4f}, p = {p_val:.4f}, dof = {dof}")

# 4. Korelasi Pearson
x = np.random.randn(100)
y = 2 * x + np.random.randn(100) * 0.5
r, p = stats.pearsonr(x, y)
print(f"\n=== Korelasi Pearson ===")
print(f"r = {r:.4f}, p = {p:.6f}")

7. Integration & ODE

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt

# Numerical integration: ∫₀¹ x² dx = 1/3
def f(x):
    return x**2

hasil, error = integrate.quad(f, 0, 1)
print(f"∫₀¹ x² dx = {hasil:.6f} (exact: {1/3:.6f}), error: {error:.2e}")

# Double integral: ∫₀¹ ∫₀¹ xy dx dy = 1/4
def f2(y, x):  # Note: inner integral variable first
    return x * y

hasil2, _ = integrate.dblquad(f2, 0, 1, 0, 1)
print(f"∫₀¹∫₀¹ xy dx dy = {hasil2:.6f} (exact: 0.250000)")

# ODE Solver: Model SIR (Epidemiology)
def sir_model(t, y, beta=0.3, gamma=0.1):
    S, I, R = y
    dSdt = -beta * S * I
    dIdt = beta * S * I - gamma * I
    dRdt = gamma * I
    return [dSdt, dIdt, dRdt]

# Kondisi awal: 99% susceptible, 1% infected, 0% recovered
y0 = [0.99, 0.01, 0]
t_span = (0, 100)
t_eval = np.linspace(0, 100, 1000)

sol = integrate.solve_ivp(sir_model, t_span, y0, t_eval=t_eval)

plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], 'b-', label='Susceptible')
plt.plot(sol.t, sol.y[1], 'r-', label='Infected')
plt.plot(sol.t, sol.y[2], 'g-', label='Recovered')
plt.xlabel('Hari')
plt.ylabel('Proporsi Populasi')
plt.title('Model SIR - Simulasi Epidemi')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('sir_model.png', dpi=150)
plt.show()

print(f"Peak Infected: {sol.y[1].max():.2%} pada hari {sol.t[sol.y[1].argmax()]:.0f}")

8. Sparse Matrices

import numpy as np
from scipy import sparse

# Matrix besar tapi mayoritas nol
n = 10000
density = 0.001  # Hanya 0.1% non-zero

# Buat sparse random matrix
dense = np.random.choice([0.0, 1.0], size=(n, n), p=[1-density, density])
sparse_matrix = sparse.csr_matrix(dense)

print(f"Dense: {dense.nbytes / 1e6:.1f} MB")
print(f"Sparse: {sparse_matrix.data.nbytes / 1e6:.1f} MB")
print(f"Ratio: {dense.nbytes / sparse_matrix.data.nbytes:.0f}x lebih hemat")

# Operasi sparse
A = sparse.random(1000, 1000, density=0.01, format='csr')
b = np.random.randn(1000)

# Sparse linear solve (jauh lebih cepat untuk matrix besar)
from scipy.sparse.linalg import spsolve, eigsh

# Solve Ax = b
x = spsolve(A, b)
print(f"Residual norm: {np.linalg.norm(A @ x - b):.2e}")

# Eigenvalues terbesar
eigenvalues, eigenvectors = eigsh(A, k=5)
print(f"Top 5 eigenvalues: {eigenvalues}")

9. Spatial & Image Processing

import numpy as np
from scipy import spatial, ndimage
import matplotlib.pyplot as plt

# KDTree untuk nearest neighbor search
np.random.seed(42)
points = np.random.randn(10000, 2)
tree = spatial.KDTree(points)

# Cari 5 titik terdekat dari query point
query = np.array([0.5, 0.5])
distances, indices = tree.query(query, k=5)

print(f"Query point: {query}")
print(f"5 nearest neighbors:")
for i, (d, idx) in enumerate(zip(distances, indices)):
    print(f"  {i+1}. Point {idx}: {points[idx]}, distance={d:.4f}")

# Convex Hull
hull = spatial.ConvexHull(points[:100])
print(f"\nConvex Hull: {hull.vertices.shape[0]} vertices")
print(f"Area: {hull.area:.4f}")

# Image processing dengan ndimage
# Buat "image" random
img = np.random.rand(100, 100)
img[30:70, 30:70] += 2  # Tambah "objek"

# Gaussian blur
img_blur = ndimage.gaussian_filter(img, sigma=3)

# Edge detection
img_sobel = ndimage.sobel(img)

# Label connected components
threshold = img > 1.5
labeled, num_features = ndimage.label(threshold)
print(f"\nDitemukan {num_features} objek terpisah")

fig, axes = plt.subplots(1, 4, figsize=(16, 4))
axes[0].imshow(img, cmap='gray')
axes[0].set_title('Original')
axes[1].imshow(img_blur, cmap='gray')
axes[1].set_title('Gaussian Blur')
axes[2].imshow(img_sobel, cmap='gray')
axes[2].set_title('Sobel Edge')
axes[3].imshow(labeled, cmap='nipy_spectral')
axes[3].set_title(f'Labeling ({num_features} objek)')
for ax in axes:
    ax.axis('off')
plt.tight_layout()
plt.savefig('ndimage.png', dpi=150)
plt.show()

10. Studi Kasus Praktis

Studi Kasus: Analisis Sinyal EKG

import numpy as np
from scipy import signal, stats
import matplotlib.pyplot as plt

# Simulasi sinyal EKG sederhana
fs = 500  # Sample rate
t = np.arange(0, 5, 1/fs)

# Sinyal EKG sintetik (QRS complex sederhana)
def ecg_wave(t, freq=1.2):
    """Generate synthetic ECG-like signal"""
    phase = (t * freq) % 1
    # P wave, QRS, T wave (simplified)
    p = 0.15 * np.exp(-((phase - 0.1)**2) / 0.002)
    qrs = 1.0 * np.exp(-((phase - 0.3)**2) / 0.0005)
    t_wave = 0.3 * np.exp(-((phase - 0.55)**2) / 0.008)
    return p + qrs + t_wave

# Sinyal EKG + noise
heart_rate = 72  # BPM
ecg_clean = ecg_wave(t, freq=heart_rate/60)
noise = 0.1 * np.random.randn(len(t))
ecg_noisy = ecg_clean + noise

# Filter: bandpass 1-40 Hz
b, a = signal.butter(4, [1, 40], btype='bandpass', fs=fs)
ecg_filtered = signal.filtfilt(b, a, ecg_noisy)

# Deteksi R-peaks (heartbeat detection)
r_peaks, _ = signal.find_peaks(ecg_filtered,
    height=0.5,
    distance=int(fs * 0.5)  # Min 0.5s between beats
)

# Hitung heart rate
rr_intervals = np.diff(r_peaks) / fs  # dalam detik
heart_rates = 60 / rr_intervals  # BPM

print(f"R-Peaks terdeteksi: {len(r_peaks)}")
print(f"Rata-rata Heart Rate: {heart_rates.mean():.1f} BPM")
print(f"HRV (Std RR): {rr_intervals.std()*1000:.1f} ms")
print(f"Min HR: {heart_rates.min():.1f}, Max HR: {heart_rates.max():.1f}")

# Visualisasi
fig, axes = plt.subplots(3, 1, figsize=(14, 8))
axes[0].plot(t[:1000], ecg_noisy[:1000], 'gray', alpha=0.7)
axes[0].set_title('Sinyal EKG Noisy')

axes[1].plot(t[:1000], ecg_filtered[:1000], 'b-')
r_in_range = r_peaks[r_peaks < 1000]
axes[1].plot(t[r_in_range], ecg_filtered[r_in_range], 'rv', markersize=10)
axes[1].set_title('EKG Filtered + R-Peaks')

axes[2].bar(range(len(heart_rates)), heart_rates, color='green', alpha=0.7)
axes[2].axhline(heart_rates.mean(), color='red', linestyle='--', label=f'Mean: {heart_rates.mean():.1f}')
axes[2].set_title('Heart Rate per Beat')
axes[2].set_xlabel('Beat Number')
axes[2].set_ylabel('BPM')
axes[2].legend()

plt.tight_layout()
plt.savefig('ecg_analysis.png', dpi=150)
plt.show()

11. Quiz Pemahaman