xbe.at

Statistics and maths / Probability and distributions

Gaussian Filters for Smoothing Time Series Data

6 min read, 14 sections

Contents
  1. Understanding Gaussian Filters for Time Series
  2. Implementing Basic Time Series Smoothing
  3. Edge Effects and Padding Strategies
  4. Adaptive Gaussian Filtering
  5. Real-time Gaussian Filtering
  6. Performance Optimization with NumPy Vectorization
  7. Multi-dimensional Gaussian Filtering
  8. Handling Missing Data in Time Series
  9. Mathematical Foundations
  10. Real-world Application - Financial Data
  11. Real-world Application - Signal Processing
  12. Adaptive Window Size Selection
  13. Performance Benchmarking
  14. Additional Resources

Understanding Gaussian Filters for Time Series

The Gaussian filter is a fundamental smoothing technique that applies a weighted average using the normal distribution. It's particularly effective for reducing noise while preserving underlying trends in time series data, making it invaluable for signal processing and data analysis.

import numpy as np
import matplotlib.pyplot as plt

def gaussian_kernel(size, sigma):
    """Generate 1D Gaussian kernel"""
    x = np.linspace(-size/2, size/2, size)
    kernel = np.exp(-x**2 / (2*sigma**2))
    return kernel / kernel.sum()

# Example kernel visualization
kernel = gaussian_kernel(25, 3)
plt.plot(kernel)
plt.title('Gaussian Kernel (σ=3)')
plt.show()

Implementing Basic Time Series Smoothing

A practical implementation of Gaussian smoothing involves convolving the input signal with a Gaussian kernel. This process effectively weights neighboring points according to their distance from the central point being smoothed.

def smooth_timeseries(data, window_size, sigma):
    kernel = gaussian_kernel(window_size, sigma)
    # Reflect padding to handle edges
    padded = np.pad(data, window_size//2, mode='reflect')
    smoothed = np.convolve(padded, kernel, mode='valid')
    return smoothed

# Generate sample noisy data
t = np.linspace(0, 10, 1000)
signal = np.sin(t) + np.random.normal(0, 0.2, len(t))
smoothed = smooth_timeseries(signal, 25, 3)

Edge Effects and Padding Strategies

Edge effects in Gaussian filtering can significantly impact results near the boundaries of the time series. Different padding strategies help mitigate these effects, each with its own trade-offs for preserving signal characteristics at the edges.

def compare_padding_strategies(data, window_size, sigma):
    padding_modes = ['reflect', 'edge', 'wrap']
    results = {}
    
    for mode in padding_modes:
        padded = np.pad(data, window_size//2, mode=mode)
        kernel = gaussian_kernel(window_size, sigma)
        smoothed = np.convolve(padded, kernel, mode='valid')
        results[mode] = smoothed
    
    return results

Adaptive Gaussian Filtering

The adaptive Gaussian filter adjusts its sigma parameter based on local signal characteristics, providing better preservation of sharp transitions while still smoothing noise in stable regions.

def adaptive_gaussian_smooth(data, window_size, base_sigma, sensitivity=0.5):
    # Calculate local variance
    local_var = np.array([np.var(data[max(0, i-window_size//2):
                                     min(len(data), i+window_size//2)])
                         for i in range(len(data))])
    
    # Adjust sigma based on local variance
    adaptive_sigma = base_sigma * (1 + sensitivity * local_var/np.max(local_var))
    
    result = np.zeros_like(data)
    for i in range(len(data)):
        kernel = gaussian_kernel(window_size, adaptive_sigma[i])
        start = max(0, i - window_size//2)
        end = min(len(data), i + window_size//2 + 1)
        window = data[start:end]
        result[i] = np.sum(window * kernel[:len(window)]) / np.sum(kernel[:len(window)])
    
    return result

Real-time Gaussian Filtering

Implementation of a causal Gaussian filter suitable for real-time applications, using only past data points to compute the smoothed value for the current time step.

def realtime_gaussian_smooth(data, window_size, sigma):
    result = np.zeros_like(data)
    kernel = gaussian_kernel(window_size, sigma)[:window_size//2 + 1]
    kernel = kernel / np.sum(kernel)  # Renormalize for causal filtering
    
    for i in range(len(data)):
        start = max(0, i - window_size//2)
        window = data[start:i+1]
        current_kernel = kernel[-len(window):]
        current_kernel = current_kernel / np.sum(current_kernel)
        result[i] = np.sum(window * current_kernel)
    
    return result

Performance Optimization with NumPy Vectorization

Optimizing Gaussian filtering implementation using NumPy's vectorized operations significantly improves computational efficiency. This approach reduces processing time for large datasets while maintaining numerical accuracy.

def vectorized_gaussian_smooth(data, window_size, sigma):
    kernel = gaussian_kernel(window_size, sigma)
    # Create a matrix of shifted data for vectorized computation
    matrix = np.zeros((len(data), window_size))
    for i in range(window_size):
        shift = i - window_size//2
        matrix[:, i] = np.roll(data, shift)
    
    # Apply kernel to all points simultaneously
    smoothed = np.sum(matrix * kernel, axis=1)
    return smoothed

Multi-dimensional Gaussian Filtering

Extending the Gaussian filter to handle multi-dimensional time series data, useful for processing multiple synchronized signals or spatiotemporal data.

def multidim_gaussian_smooth(data, window_size, sigma):
    """
    data: shape (n_samples, n_dimensions)
    """
    kernel = gaussian_kernel(window_size, sigma)
    smoothed = np.zeros_like(data)
    
    for dim in range(data.shape[1]):
        padded = np.pad(data[:, dim], window_size//2, mode='reflect')
        smoothed[:, dim] = np.convolve(padded, kernel, mode='valid')
    
    return smoothed

Handling Missing Data in Time Series

A robust implementation of Gaussian filtering that properly handles missing values (NaN) in the time series, essential for real-world applications with incomplete data.

def robust_gaussian_smooth(data, window_size, sigma):
    kernel = gaussian_kernel(window_size, sigma)
    result = np.zeros_like(data)
    
    for i in range(len(data)):
        start = max(0, i - window_size//2)
        end = min(len(data), i + window_size//2 + 1)
        window = data[start:end]
        current_kernel = kernel[window_size//2-(i-start):window_size//2+(end-i)]
        
        # Handle NaN values
        valid_mask = ~np.isnan(window)
        if np.any(valid_mask):
            current_kernel = current_kernel * valid_mask
            current_kernel = current_kernel / np.sum(current_kernel)
            result[i] = np.sum(window[valid_mask] * current_kernel[valid_mask])
        else:
            result[i] = np.nan
            
    return result

Mathematical Foundations

The theoretical underpinnings of Gaussian filtering, including the mathematical formulation and its properties, presented with key equations and their practical implications.

# Mathematical formulation of Gaussian filter
"""
The Gaussian function in 1D:
$$G(x) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{x^2}{2\sigma^2}}$$

Discrete convolution formula:
$$y[n] = \sum_{k=-\infty}^{\infty} h[k]x[n-k]$$

Where h[k] is the Gaussian kernel:
$$h[k] = G(k\Delta t)$$
"""

Real-world Application - Financial Data

Implementation of Gaussian filtering for smoothing financial time series data, demonstrating practical application in market trend analysis.

def analyze_financial_timeseries(prices, window_size=25, sigma=3):
    # Calculate returns
    returns = np.diff(np.log(prices))
    
    # Apply Gaussian smoothing
    smoothed_returns = smooth_timeseries(returns, window_size, sigma)
    
    # Calculate volatility estimation
    volatility = np.zeros_like(returns)
    for i in range(len(returns)):
        window = returns[max(0, i-window_size):i+1]
        volatility[i] = np.std(window) * np.sqrt(252)  # Annualized
    
    return {
        'smoothed_returns': smoothed_returns,
        'volatility': volatility
    }

# Example usage with sample data
prices = np.random.lognormal(0.0001, 0.02, 1000).cumprod()
results = analyze_financial_timeseries(prices)

Real-world Application - Signal Processing

A comprehensive example of Gaussian filtering applied to signal processing, demonstrating noise reduction in sensor data while preserving important signal features.

def process_sensor_signal(signal_data, sampling_rate, cutoff_freq):
    # Calculate appropriate sigma based on cutoff frequency
    sigma = sampling_rate / (2 * np.pi * cutoff_freq)
    
    # Window size should be ~6 sigma for 99.7% of Gaussian distribution
    window_size = int(6 * sigma)
    if window_size % 2 == 0:
        window_size += 1
    
    # Apply filtering with frequency-based parameters
    filtered_signal = smooth_timeseries(signal_data, window_size, sigma)
    
    # Calculate signal-to-noise ratio improvement
    noise_before = np.std(signal_data - np.mean(signal_data))
    noise_after = np.std(filtered_signal - np.mean(filtered_signal))
    snr_improvement = 20 * np.log10(noise_before / noise_after)
    
    return filtered_signal, snr_improvement

Adaptive Window Size Selection

Implementation of an intelligent window size selection algorithm that automatically determines optimal Gaussian filter parameters based on signal characteristics.

def auto_window_size(data, target_smoothness=0.95):
    # Calculate signal properties
    signal_range = np.max(data) - np.min(data)
    noise_estimate = np.median(np.abs(np.diff(data))) * 1.4826
    
    def evaluate_smoothness(window_size):
        smoothed = smooth_timeseries(data, window_size, window_size/6)
        residuals = data - smoothed
        smoothness = 1 - (np.std(residuals) / signal_range)
        return smoothness
    
    # Binary search for optimal window size
    left, right = 3, len(data)//4
    while left < right:
        mid = (left + right) // 2
        if mid % 2 == 0:
            mid += 1
        smoothness = evaluate_smoothness(mid)
        
        if abs(smoothness - target_smoothness) < 0.01:
            return mid
        elif smoothness < target_smoothness:
            left = mid + 2
        else:
            right = mid - 2
            
    return left

Performance Benchmarking

A comprehensive benchmarking suite to evaluate different Gaussian filtering implementations in terms of speed, accuracy, and memory usage.

def benchmark_gaussian_implementations(data_size=10000, n_trials=10):
    # Generate test data
    np.random.seed(42)
    data = np.random.randn(data_size)
    
    implementations = {
        'basic': smooth_timeseries,
        'vectorized': vectorized_gaussian_smooth,
        'adaptive': adaptive_gaussian_smooth,
        'realtime': realtime_gaussian_smooth
    }
    
    results = {}
    for name, func in implementations.items():
        times = []
        memory_usage = []
        
        for _ in range(n_trials):
            start_time = time.time()
            _ = func(data, 25, 3)
            times.append(time.time() - start_time)
            
            # Measure peak memory usage
            tracemalloc.start()
            _ = func(data, 25, 3)
            current, peak = tracemalloc.get_traced_memory()
            memory_usage.append(peak / 1024)  # Convert to KB
            tracemalloc.stop()
            
        results[name] = {
            'avg_time': np.mean(times),
            'std_time': np.std(times),
            'avg_memory': np.mean(memory_usage)
        }
        
    return results

Additional Resources

More in Probability and distributions

  1. Geometric Distribution Explained with Python6 min

    The Geometric Distribution is a discrete probability distribution that models the number of trials needed to achieve the first success in a sequence of...

  2. Explaining Normal Distribution in Data Science10 min

    The Normal Distribution, also known as the Gaussian distribution, is a symmetrical, bell-shaped curve that represents the probability distribution of many...

  3. Limitations of Normal Distribution in Nature7 min

    The normal distribution, while widely used, is not always the best model for real-world phenomena. Many natural processes are asymmetric or involve...

  4. Essentials of Probability and Statistics5 min

    Probability is the branch of mathematics that deals with the likelihood of events occurring. It forms the foundation for statistical analysis and...

  5. Master Probability Distributions with Python7 min

    Probability distributions are mathematical functions that describe the likelihood of different outcomes in a random event. They are fundamental to statistics...

  6. Empirical Rule and Normal Distribution in Python9 min

    The Empirical Rule, also known as the 68-95-99.7 rule, is a fundamental concept in statistics that describes the distribution of data in a normal...