Mathematics of Industrial Signal Processing

What Is Digital Signal Processing and Why Does It Matter in Industry?

Every sensor in a factory -- temperature, pressure, vibration, flow -- produces an electrical signal carrying information. But that signal is always mixed with noise: electromagnetic interference from motors, mechanical vibrations, or power line distortion.

Digital Signal Processing (DSP) is the set of mathematical tools that separates useful information from noise, extracts hidden patterns, and compresses data for transmission and storage.

In this lesson we cover the mathematical foundations of DSP with practical examples from the factory floor.

Sampling Theorem: The Basis of Analog-to-Digital Conversion

When converting an analog signal to digital, we take samples at regular intervals. The fundamental question is: how many samples per second do we need?

The Nyquist-Shannon theorem states:

fs >= 2 * fmax

where fs is the sampling frequency and fmax is the highest frequency in the signal. The frequency fs/2 is called the Nyquist frequency.

Practical example: A vibration sensor on a ball bearing measures frequencies up to 5 kHz:

# Minimum sampling frequency
f_max = 5000  # Hz (highest vibration frequency)
fs_min = 2 * f_max  # 10,000 Hz

# In practice use at least 2.5x for safety margin
fs_practical = int(2.5 * f_max)  # 12,500 Hz

print(f"Theoretical minimum: {fs_min} Hz")
print(f"Practical sampling rate: {fs_practical} Hz")
print(f"Time between samples: {1/fs_practical * 1e6:.1f} us")

Aliasing: The Fatal Error

What happens if we sample below the Nyquist rate? High frequencies appear folded as false low frequencies -- this is called aliasing.

Imagine a car wheel in a film appearing to spin backward -- the exact same phenomenon. The camera frame rate (sampling frequency) is slower than the wheel rotation (signal frequency).

import numpy as np

# Real signal: 800 Hz
f_signal = 800  # Hz

# Sampling at 1000 Hz (less than 2 * 800)
fs_bad = 1000
t_bad = np.arange(0, 0.01, 1/fs_bad)
samples_bad = np.sin(2 * np.pi * f_signal * t_bad)

# Alias frequency caused by undersampling
f_alias = abs(f_signal - fs_bad)  # 200 Hz!
print(f"True frequency: {f_signal} Hz")
print(f"Sampling rate: {fs_bad} Hz")
print(f"Alias frequency: {f_alias} Hz")
print("The signal appears to be 200 Hz instead of 800 Hz!")

# Solution: Anti-aliasing filter before ADC, then adequate sampling
fs_good = 2000  # greater than 2 * 800
print(f"\nCorrect sampling rate: {fs_good} Hz (> 2 x {f_signal})")

Practical solution: Always use a low-pass anti-aliasing filter before the analog-to-digital converter (ADC) to remove frequencies above Nyquist.

The Z-Transform: A Tool for Analyzing Digital Filters

Just as the Laplace transform simplifies analysis of analog systems, the Z-transform does the same for digital systems. It converts difference equations into algebraic equations.

X(z) = sum of x[n] * z^(-n)    (n from 0 to infinity)

where z is a complex variable. Key properties:

Why does it matter? The transfer function of any digital filter is written in terms of z:

H(z) = (b0 + b1*z^(-1) + b2*z^(-2)) / (1 + a1*z^(-1) + a2*z^(-2))

FIR Filter: Finite Impulse Response

An FIR filter depends only on current and past input samples:

y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] + ... + bN*x[n-N]

Advantages:

• Always stable (cannot blow up)
• Can be designed with linear phase -- does not distort the signal shape
• Easy to implement on microcontrollers

Disadvantages: Requires more coefficients than IIR for the same sharpness.

import numpy as np

def fir_filter(x, coefficients):
    """Simple FIR filter"""
    N = len(coefficients)
    y = np.zeros(len(x))
    for n in range(len(x)):
        for k in range(N):
            if n - k >= 0:
                y[n] += coefficients[k] * x[n - k]
    return y

# Moving average filter (simplest FIR) -- 5 taps
coefficients = [0.2, 0.2, 0.2, 0.2, 0.2]

# Temperature sensor signal with noise
np.random.seed(42)
t = np.arange(100)
signal_clean = 75.0 + 0.1 * t  # gradually rising temperature
noise = np.random.normal(0, 2, 100)  # random noise
signal_noisy = signal_clean + noise

# Apply the filter
signal_filtered = fir_filter(signal_noisy, coefficients)

print(f"Raw reading (sample 50): {signal_noisy[50]:.2f} C")
print(f"After filtering: {signal_filtered[50]:.2f} C")
print(f"True value: {signal_clean[50]:.2f} C")

IIR Filter: Infinite Impulse Response

An IIR filter depends on both input samples and previous output samples:

y[n] = b0*x[n] + b1*x[n-1] - a1*y[n-1] - a2*y[n-2]

Advantages: Fewer coefficients for the same performance -- ideal for resource-constrained systems.

Disadvantages: Can become unstable if not designed carefully.

def iir_lowpass(x, alpha=0.1):
    """First-order IIR low-pass filter
    alpha: smoothing factor (0 < alpha < 1)
    smaller alpha = stronger smoothing but slower response
    """
    y = np.zeros(len(x))
    y[0] = x[0]
    for n in range(1, len(x)):
        y[n] = alpha * x[n] + (1 - alpha) * y[n-1]
    return y

# Compare different alpha values
for alpha in [0.05, 0.1, 0.3]:
    filtered = iir_lowpass(signal_noisy, alpha)
    error = np.mean(np.abs(filtered[10:] - signal_clean[10:]))
    print(f"alpha={alpha:.2f}: mean error = {error:.3f} C")

Typical output:

alpha=0.05: mean error = 1.102 C
alpha=0.10: mean error = 0.873 C
alpha=0.30: mean error = 1.241 C

The best alpha depends on the signal and noise characteristics -- practical experimentation is essential.

Windowing: Improving Spectral Analysis

When analyzing a signal using the Fast Fourier Transform (FFT), we cut the signal into finite segments. This truncation causes spectral leakage -- distortion in the frequency spectrum.

Solution: Multiply the signal by a window function before FFT:

import numpy as np

N = 256  # number of samples
n = np.arange(N)

# Test signal: two close frequencies
fs = 1000
f1, f2 = 100, 110  # only 10 Hz apart
signal = np.sin(2*np.pi*f1/fs*n) + 0.5*np.sin(2*np.pi*f2/fs*n)

# Without window (rectangular)
fft_rect = np.abs(np.fft.fft(signal))[:N//2]

# With Hamming window
window = np.hamming(N)
fft_hamming = np.abs(np.fft.fft(signal * window))[:N//2]

freqs = np.fft.fftfreq(N, 1/fs)[:N//2]

# Compare peaks
peak_rect = freqs[np.argsort(fft_rect)[-2:]]
peak_hamming = freqs[np.argsort(fft_hamming)[-2:]]

print(f"Peaks without window: {sorted(peak_rect)} Hz")
print(f"Peaks with Hamming: {sorted(peak_hamming)} Hz")
print(f"True frequencies: [{f1}, {f2}] Hz")

Industrial Application: Bearing Fault Detection via Vibration Analysis

In a factory, a vibration sensor on a motor measures a complex signal. Using DSP tools we can detect ball bearing faults early:

import numpy as np

# Simulate vibration signal from a ball bearing
np.random.seed(42)
fs = 10000  # sampling rate 10 kHz
t = np.arange(0, 1, 1/fs)  # one second

# Normal vibration (50 Hz from motor rotation)
vibration = 2.0 * np.sin(2*np.pi*50*t)

# Add outer race fault frequency (BPFO = 162 Hz)
# This frequency indicates early outer race damage
fault_signal = 0.3 * np.sin(2*np.pi*162*t)
vibration += fault_signal

# Noise
vibration += np.random.normal(0, 0.5, len(t))

# FFT analysis with Hamming window
N = len(vibration)
window = np.hamming(N)
fft_result = np.abs(np.fft.fft(vibration * window))[:N//2]
freqs = np.fft.fftfreq(N, 1/fs)[:N//2]

# Find peaks above threshold
threshold = np.mean(fft_result) + 3 * np.std(fft_result)
peaks_idx = np.where(fft_result > threshold)[0]
peak_freqs = freqs[peaks_idx]

print("Detected frequencies (above threshold):")
for f in sorted(set(np.round(peak_freqs, 0))):
    if f > 0:
        print(f"  {f:.0f} Hz", end="")
        if abs(f - 50) < 5:
            print(" (motor rotation - normal)")
        elif abs(f - 162) < 5:
            print(" (BPFO - WARNING: bearing fault!)")
        else:
            print()

# Severity indicator
rms_total = np.sqrt(np.mean(vibration**2))
print(f"\nOverall RMS: {rms_total:.3f} mm/s")
if rms_total < 2.8:
    print("Condition: Good")
elif rms_total < 7.1:
    print("Condition: Acceptable - monitoring required")
else:
    print("Condition: Dangerous - immediate maintenance!")

Summary and Comparison

Practical tip: On the factory floor, always start with a simple moving average (FIR). If that is not sufficient, step up to a first-order IIR (Exponential Moving Average). Do not overcomplicate things unless the application demands it -- the simplest filter that meets requirements is the best choice.