### Railway Accelerometer Analysis with VKF

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

A comprehensive example simulating and analyzing railway axle box accelerometer data using VKF to extract wheel and track vibration components. This demonstrates VKF's application in real-world scenarios with time-varying frequencies.

```python
import numpy as np
import matplotlib.pyplot as plt
from VoldKalmanFilter import vkf

# Simulate railway accelerometer signal
fs = 12000      # 12 kHz sampling (typical for ABA)
T = 60          # 60 seconds of data
dt = 1/fs
t = np.arange(0, T + dt, dt)
N = t.size

# Train speed profile (accelerating from 20 to 40 m/s)
speed = 20 + 20 * t/T  # m/s

# Wheel diameter and circumference
wheel_diameter = 0.92  # meters (typical passenger car wheel)
wheel_circumference = np.pi * wheel_diameter

# Wheel rotation frequency (speed / circumference)
f_wheel = (speed / wheel_circumference).reshape(-1, 1)

# Wheel harmonics (out-of-roundness creates harmonics)
f_wheel_1x = f_wheel           # 1st harmonic
f_wheel_2x = 2 * f_wheel       # 2nd harmonic
f_wheel_3x = 3 * f_wheel       # 3rd harmonic

# Simulate wheel vibration components
A_1x = 2.0   # Amplitude for 1st harmonic (wheel flat)
A_2x = 0.8   # Amplitude for 2nd harmonic (ovalization)
A_3x = 0.3   # Amplitude for 3rd harmonic

y_1x = A_1x * np.cos(2*np.pi * np.cumsum(f_wheel_1x) * dt).flatten()
y_2x = A_2x * np.cos(2*np.pi * np.cumsum(f_wheel_2x) * dt).flatten()
y_3x = A_3x * np.cos(2*np.pi * np.cumsum(f_wheel_3x) * dt).flatten()

# Add track excitation and noise
track_noise = 0.5 * np.random.randn(N)
measurement_noise = 0.2 * np.random.randn(N)

# Combined accelerometer signal
y_aba = y_1x + y_2x + y_3x + track_noise + measurement_noise
```

--------------------------------

### VKF with Boundary Conditions

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Apply advanced filter orders with boundary conditions by specifying an array for the `p` parameter. This forces the envelope amplitude and its derivatives to be zero at the signal boundaries, useful for avoiding edge effects.

```python
import numpy as np
from VoldKalmanFilter import vkf

fs = 8000
T = 10
dt = 1/fs
t = np.arange(0, T + dt, dt)
N = t.size

# Create a burst signal (transient)
f0 = 200 * np.ones((N, 1))
envelope = np.exp(-((t - T/2)**2) / 2)  # Gaussian envelope
y_signal = envelope * np.cos(2*np.pi * np.cumsum(f0) * dt).flatten()
noise = 0.1 * np.random.randn(N)
y = y_signal + noise

# Apply 2nd order filtering with boundary conditions
# p = [2, 0, 1] means:
#   - 2nd order main filter
#   - Force amplitude to zero at boundaries (order 0)
#   - Force first derivative to zero at boundaries (order 1)
(x, c, r) = vkf(y, fs, f0, p=[2, 0, 1], bw=1.0)

xreal = np.real(x * c)
extracted_envelope = np.abs(x).flatten()

# Verify boundary conditions
print(f"Start amplitude: {extracted_envelope[0]:.6f}")
print(f"End amplitude: {extracted_envelope[-1]:.6f}")
print(f"Peak amplitude: {np.max(extracted_envelope):.4f}")
```

--------------------------------

### Apply Vold-Kalman Filter to Synthetic Signal

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Demonstrates applying the vkf function to a synthetic signal with multiple stationary and non-stationary components. Requires numpy and the vkf function. The output includes the complex envelope, phasor, and selectivity arrays, which can be used to reconstruct filtered signals and analyze component properties.

```python
import numpy as np
from VoldKalmanFilter import vkf

# Setup signal parameters
fs = 12000  # Sampling frequency (Hz)
T = 50      # Duration (seconds)
dt = 1/fs
t = np.arange(0, T + dt, dt)
N = t.size

# Create synthetic signal with three components
# Component 1: Non-stationary frequency with varying amplitude
A0 = 2 * np.linspace(0.5, 1, N)
f0 = (fs/16 - fs/100 * np.cos(2*np.pi*1.5*t/T)).reshape(-1, 1)
phi0 = 2*np.pi * np.cumsum(f0) * dt
y0 = A0 * np.cos(phi0.flatten())

# Component 2: Another non-stationary component
A1 = 0.7 + 0.4 * np.cos(t*0.04*np.pi*2) * np.sin(t*0.02*np.pi*2)
f1 = (fs/50 - fs/100 * np.cos(2*np.pi*t/T/2)).reshape(-1, 1)
phi1 = 2*np.pi * np.cumsum(f1) * dt
y1 = A1 * np.cos(phi1.flatten())

# Component 3: Stationary component at 500 Hz
As1 = np.ones(N)
fs1 = 500 * np.ones((N, 1))
phis1 = 2*np.pi * np.cumsum(fs1) * dt
ys1 = As1 * np.cos(phis1.flatten())

# Add white noise
noise = (np.random.rand(N) - 0.5) * 2.6

# Mixed signal
y = y0 + y1 + ys1 + noise

# Apply Vold-Kalman Filter
# Parameters:
#   y: input signal
#   fs: sampling frequency
#   f: frequency array (N x K) for K orders
#   p: filter order (1-4, higher = sharper roll-off at -40dB/decade per order)
#   bw: bandwidth in Hz (scalar for constant, array for time-varying)
p = 2           # 2nd order filter
bw = fs/8000    # Bandwidth of 1.5 Hz

(x, c, r) = vkf(y, fs, np.hstack([f0, f1, fs1]), p=p, bw=bw)
# Returns:
#   x: complex envelope array (N x K)
#   c: phasor array (N x K)
#   r: selectivity array

# Reconstruct filtered real signals
xreal = np.real(x * c)

# Extract individual components
component_1 = xreal[:, 0]  # First frequency component
component_2 = xreal[:, 1]  # Second frequency component
component_3 = xreal[:, 2]  # Third frequency component (500 Hz stationary)

# Get amplitude and phase of each component
amplitude = np.abs(x)      # Instantaneous amplitude envelope
phase = np.angle(x)        # Instantaneous phase

# Reveal residual noise
residual = y - np.sum(xreal, axis=1)

# Calculate Mean Squared Error
mse = np.mean((y0 + y1 + ys1 - np.sum(xreal, axis=1))**2)
print(f"MSE: {mse:.6f}")
```

--------------------------------

### Sparse Matrix Row Repetition Utility

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Use `lil_repeat` to repeat rows in LIL sparse matrices. This utility is used internally for constructing filter matrices and can be helpful for custom sparse matrix operations.

```python
import numpy as np
from scipy import sparse
from VoldKalmanFilter import lil_repeat

# Create a small sparse matrix in LIL format
data = np.array([[1, 0, 2], [0, 3, 0], [4, 0, 5]])
S = sparse.lil_matrix(data)

print("Original matrix:")
print(S.toarray())
# [[1 0 2]
#  [0 3 0]
#  [4 0 5]]

# Repeat each row 2 times
repeated = lil_repeat(S, 2)
print("\nRepeated matrix (each row 2x):")
print(repeated.toarray())
# [[1 0 2]
#  [1 0 2]
#  [0 3 0]
#  [0 3 0]
#  [4 0 5]
#  [4 0 5]]

# Variable repetition per row
S2 = sparse.lil_matrix(np.array([[1, 2], [3, 4], [5, 6]]))
repeated_var = lil_repeat(S2, [1, 3, 2])  # Row 0: 1x, Row 1: 3x, Row 2: 2x
print("\nVariable repetition:")
print(repeated_var.toarray())
```

--------------------------------

### VKF with Time-Varying Bandwidth

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Use a time-varying bandwidth for adaptive filtering that can track changing signal characteristics. Pass an array to the `bw` parameter instead of a scalar.

```python
import numpy as np
from VoldKalmanFilter import vkf

fs = 12000
T = 30
dt = 1/fs
t = np.arange(0, T + dt, dt)
N = t.size

# Signal with frequency that changes over time (chirp)
f_inst = (100 + 50 * t/T).reshape(-1, 1)  # 100 Hz to 150 Hz sweep
y_signal = np.cos(2*np.pi * np.cumsum(f_inst) * dt).flatten()
noise = 0.3 * np.random.randn(N)
y = y_signal + noise

# Time-varying bandwidth: narrower at start, wider at end
# This can help track rapidly changing frequencies
bw_varying = np.linspace(0.5, 3.0, N).reshape(-1, 1)

(x, c, r) = vkf(y, fs, f_inst, p=2, bw=bw_varying)

xreal = np.real(x * c)
amplitude_envelope = np.abs(x)

print(f"Mean amplitude: {np.mean(amplitude_envelope):.4f}")
print(f"Amplitude std: {np.std(amplitude_envelope):.4f}")
```

--------------------------------

### Print Wheel Condition Analysis Results

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Prints the calculated RMS values for the 1x, 2x, and 3x wheel harmonics, along with the RMS value for the residual signal representing track excitation. Formats output for clarity.

```python
print("Wheel Condition Analysis Results:")
print(f"1x Harmonic RMS: {rms_1x:.4f} (wheel flat indicator)")
print(f"2x Harmonic RMS: {rms_2x:.4f} (ovalization indicator)")
print(f"3x Harmonic RMS: {rms_3x:.4f} (polygonization indicator)")
```

--------------------------------

### Reconstruct Filtered Signals from VKF Output

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Reconstructs the filtered signals from the complex envelopes (x) and phase information (c) obtained from the vkf function. Extracts the real part to obtain the filtered time-domain signals for each harmonic.

```python
# Reconstruct filtered signals
xreal = np.real(x * c)
wheel_1x_filtered = xreal[:, 0]
wheel_2x_filtered = xreal[:, 1]
wheel_3x_filtered = xreal[:, 2]
```

--------------------------------

### Calculate RMS Values for Wheel Condition Analysis

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Computes the Root Mean Square (RMS) values of the extracted amplitude envelopes. These RMS values provide quantitative indicators for wheel condition, such as flat spots, ovalization, and polygonization.

```python
# Calculate RMS values as condition indicators
rms_1x = np.sqrt(np.mean(amplitude_1x**2))
rms_2x = np.sqrt(np.mean(amplitude_2x**2))
rms_3x = np.sqrt(np.mean(amplitude_3x**2))
```

--------------------------------

### Apply VKF for Wheel Harmonic Extraction

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Applies the Variational Kalman Filter (VKF) to extract specific harmonic frequencies from a signal. Configure filter order (p) and bandwidth (bw) for precise tracking. Requires NumPy for frequency matrix construction.

```python
frequency_matrix = np.hstack([f_wheel_1x, f_wheel_2x, f_wheel_3x])

(x, c, r) = vkf(
    y_aba,
    fs,
    frequency_matrix,
    p=2,        # 2nd order filter
    bw=0.5      # Narrow bandwidth for precise tracking
)
```

--------------------------------

### Apply Single-Order Vold-Kalman Filter for Large Datasets

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Utilizes the single-order algorithm (`multiorder=0`) of the vkf function for improved computational performance on high sample rate or long time series data. This approach processes each order independently, offering a trade-off between accuracy and speed. Requires numpy and the vkf function.

```python
import numpy as np
from VoldKalmanFilter import vkf

# High sample rate signal
fs = 48000  # 48 kHz sampling
T = 120     # 2 minutes of data
dt = 1/fs
t = np.arange(0, T + dt, dt)
N = t.size

# Create signal with known frequency components
f1 = 100 * np.ones((N, 1))   # 100 Hz component
f2 = 250 * np.ones((N, 1))   # 250 Hz component
f3 = 500 * np.ones((N, 1))   # 500 Hz component

A1, A2, A3 = 1.0, 0.5, 0.3
y1 = A1 * np.cos(2*np.pi * np.cumsum(f1) * dt).flatten()
y2 = A2 * np.cos(2*np.pi * np.cumsum(f2) * dt).flatten()
y3 = A3 * np.cos(2*np.pi * np.cumsum(f3) * dt).flatten()
noise = 0.1 * np.random.randn(N)

y = y1 + y2 + y3 + noise

# Use single-order algorithm for large datasets
(x, c, r) = vkf(
    y,
    fs,
    np.hstack([f1, f2, f3]),
    p=2,              # 2nd order filter
    bw=2.0,           # 2 Hz bandwidth
    multiorder=0      # Single-order algorithm (faster for large data)
)

xreal = np.real(x * c)
print(f"Extracted amplitudes: {np.mean(np.abs(x), axis=0)}")
# Expected output close to: [1.0, 0.5, 0.3]
```

--------------------------------

### Extract Amplitude Envelopes for Condition Indicators

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Calculates the amplitude envelopes of the filtered signals by taking the absolute value of the complex envelopes. These amplitudes serve as indicators of wheel condition.

```python
# Extract amplitude envelopes (wheel condition indicators)
amplitude_1x = np.abs(x[:, 0])
amplitude_2x = np.abs(x[:, 1])
amplitude_3x = np.abs(x[:, 2])
```

--------------------------------

### Calculate Track Excitation Residual and RMS

Source: https://context7.com/cyprienhoelzl/pyvkf/llms.txt

Calculates the residual signal by subtracting the sum of the reconstructed filtered signals (xreal) from the original input signal (y_aba). Computes the RMS of this residual to quantify track excitation.

```python
# Residual contains track excitation
residual = y_aba - np.sum(xreal, axis=1)
track_rms = np.sqrt(np.mean(residual**2))
print(f"Track excitation RMS: {track_rms:.4f}")
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.