Signal Processing & Acoustic Data Analysis: Practical Python Guide

Signal processing is one of the most fundamental disciplines in engineering. It spans a wide range of applications — from audio signals to communication systems, medical devices to underwater acoustics. In this article, we will explore Fourier transforms, spectrogram analysis, noise filtering techniques, and their practical Python implementations. I will also share insights from my master's thesis on underwater acoustic message reconstruction.

Signal Processing Fundamentals

What Is a Signal?

A signal is a mathematical representation of a quantity that varies with time or space. Signals fall into two main categories:

• Analog (Continuous) Signals: Defined on a continuous time axis (sound waves, radio signals)
• Digital (Discrete) Signals: Sampled at specific time points (digital audio, sensor data)

Sampling and the Nyquist Theorem

When converting an analog signal to digital, the sampling frequency (fs) must be at least twice the highest frequency (fmax) in the signal:

*fs >= 2 fmax**

Violating this rule causes aliasing, which irreversibly distorts the signal.

Fourier Transform and FFT

Fast Fourier Transform (FFT)

The FFT converts a time-domain signal into the frequency domain, revealing which frequency components make up the signal:

import numpy as np
import matplotlib.pyplot as plt

# Sampling parameters
fs = 8000  # Sampling frequency (Hz)
T = 1.0    # Duration (seconds)
N = int(fs * T)  # Total number of samples
t = np.linspace(0, T, N, endpoint=False)

# Create composite signal: 440 Hz + 1000 Hz + noise
f1, f2 = 440, 1000  # Frequencies (Hz)
signal = 0.7 * np.sin(2 * np.pi * f1 * t) + 0.3 * np.sin(2 * np.pi * f2 * t)
noise = 0.2 * np.random.randn(N)
noisy_signal = signal + noise

# Compute FFT
fft_result = np.fft.rfft(noisy_signal)
frequencies = np.fft.rfftfreq(N, d=1/fs)
magnitude = 2.0 / N * np.abs(fft_result)

# Plot frequency spectrum
plt.figure(figsize=(12, 4))
plt.plot(frequencies, magnitude)
plt.title('Frequency Spectrum (FFT)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim(0, 2000)
plt.grid(True, alpha=0.3)
plt.axvline(x=440, color='r', linestyle='--', alpha=0.7, label='440 Hz')
plt.axvline(x=1000, color='g', linestyle='--', alpha=0.7, label='1000 Hz')
plt.legend()
plt.tight_layout()
plt.savefig('fft_spectrum.png', dpi=150)
plt.show()

Spectrogram Analysis

A spectrogram visualizes how the frequency content of a signal changes over time. It is generated using the Short-Time Fourier Transform (STFT):

from scipy import signal as scipy_signal

# Create a chirp signal (time-varying frequency)
fs = 16000
t = np.linspace(0, 2, 2 * fs)
# Signal sweeping from 100 Hz to 4000 Hz
chirp = scipy_signal.chirp(t, f0=100, f1=4000, t1=2, method='linear')

# Compute and plot spectrogram
fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# Time domain
axes[0].plot(t, chirp, linewidth=0.5)
axes[0].set_title('Time-Domain Signal')
axes[0].set_xlabel('Time (s)')
axes[0].set_ylabel('Amplitude')
axes[0].set_xlim(0, 2)

# Spectrogram
f, t_spec, Sxx = scipy_signal.spectrogram(chirp, fs=fs, nperseg=512, noverlap=384)
axes[1].pcolormesh(t_spec, f, 10 * np.log10(Sxx + 1e-10), shading='gouraud', cmap='viridis')
axes[1].set_title('Spectrogram')
axes[1].set_xlabel('Time (s)')
axes[1].set_ylabel('Frequency (Hz)')
axes[1].set_ylim(0, 5000)
plt.colorbar(axes[1].collections[0], ax=axes[1], label='Power (dB)')

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

Digital Filtering

Noise Removal

Real-world signals typically contain noise. Digital filters can remove unwanted frequency components:

from scipy.signal import butter, filtfilt, iirnotch

def bandpass_filter(signal_data, low_cut, high_cut, fs, order=4):
    """Apply a Butterworth bandpass filter."""
    nyq = 0.5 * fs
    low = low_cut / nyq
    high = high_cut / nyq
    b, a = butter(order, [low, high], btype='band')
    filtered = filtfilt(b, a, signal_data)
    return filtered

def notch_filter(signal_data, target_freq, fs, quality_factor=30):
    """Apply a notch filter to suppress a specific frequency."""
    b, a = iirnotch(target_freq, quality_factor, fs)
    filtered = filtfilt(b, a, signal_data)
    return filtered

# Filter applications
# 300-3400 Hz bandpass (for voice communication)
voice_filtered = bandpass_filter(noisy_signal, 300, 3400, fs)

# Suppress 50 Hz mains hum
mains_clean = notch_filter(noisy_signal, 50, fs)

# Compare results
fig, axes = plt.subplots(3, 1, figsize=(12, 8))
axes[0].plot(t[:1000], noisy_signal[:1000])
axes[0].set_title('Original (Noisy) Signal')
axes[1].plot(t[:1000], voice_filtered[:1000])
axes[1].set_title('After Bandpass Filter')
axes[2].plot(t[:1000], mains_clean[:1000])
axes[2].set_title('After Notch Filter')
for ax in axes:
    ax.set_xlabel('Time (s)')
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('filtering_results.png', dpi=150)
plt.show()

Wavelet Transform

While the Fourier transform reveals the frequency content of a signal, it loses time information. STFT partially addresses this but makes a fixed tradeoff between time and frequency resolution due to its constant window size. The wavelet transform overcomes this limitation: it provides good frequency resolution at low frequencies and good time resolution at high frequencies.

This property makes the wavelet transform particularly powerful for analyzing signals containing transient events. Earthquake signals, speech, music, and biomedical signals whose frequency content changes over time are ideal candidates for wavelet analysis.

import pywt

# Continuous Wavelet Transform (CWT)
scales = np.arange(1, 128)
coefficients, frequencies_wt = pywt.cwt(noisy_signal[:2000], scales, 'morl', sampling_period=1/fs)

plt.figure(figsize=(12, 6))
plt.imshow(np.abs(coefficients), aspect='auto', cmap='jet',
           extent=[0, 2000/fs, frequencies_wt[-1], frequencies_wt[0]])
plt.title('Continuous Wavelet Transform (Morlet)')
plt.xlabel('Time (s)')
plt.ylabel('Frequency (Hz)')
plt.colorbar(label='Magnitude')
plt.tight_layout()
plt.savefig('wavelet_transform.png', dpi=150)
plt.show()

# Noise reduction with Discrete Wavelet Transform
def wavelet_denoise(signal_data, wavelet='db4', level=4, threshold_multiplier=1.0):
    """Noise reduction using wavelet thresholding."""
    coeffs = pywt.wavedec(signal_data, wavelet, level=level)

    # Apply soft thresholding to detail coefficients
    sigma = np.median(np.abs(coeffs[-1])) / 0.6745
    threshold = threshold_multiplier * sigma * np.sqrt(2 * np.log(len(signal_data)))

    clean_coeffs = [coeffs[0]]  # Keep approximation coefficients
    for detail in coeffs[1:]:
        clean_coeffs.append(pywt.threshold(detail, threshold, mode='soft'))

    clean_signal = pywt.waverec(clean_coeffs, wavelet)
    return clean_signal

denoised = wavelet_denoise(noisy_signal)

Mel-Frequency Cepstral Coefficients (MFCC)

MFCC is the most widely used feature extraction method in audio and speech processing. It mimics the human auditory system's frequency perception: high resolution at low frequencies and low resolution at high frequencies. This characteristic has made MFCC the standard feature vector for tasks like speech recognition, speaker verification, and audio classification.

MFCC computation involves several steps: signal framing, window function application, FFT computation, Mel filter bank filtering, logarithmic power computation, and finally the discrete cosine transform (DCT).

import librosa
import librosa.display

# Compute MFCC for audio signal
# Create sample signal (in real applications, loaded from audio file)
sr = 16000  # Sample rate
duration = 2.0
t_audio = np.linspace(0, duration, int(sr * duration), endpoint=False)
audio_signal = 0.5 * np.sin(2 * np.pi * 440 * t_audio) * np.exp(-t_audio)

# MFCC extraction
mfcc_coefficients = librosa.feature.mfcc(
    y=audio_signal, sr=sr,
    n_mfcc=13,           # 13 MFCC coefficients (standard)
    n_fft=2048,           # FFT window size
    hop_length=512,       # Window hop size
    n_mels=40             # Number of Mel filters
)

# Delta and delta-delta coefficients (temporal change information)
mfcc_delta = librosa.feature.delta(mfcc_coefficients)
mfcc_delta2 = librosa.feature.delta(mfcc_coefficients, order=2)

# Concatenate all features
features = np.concatenate([mfcc_coefficients, mfcc_delta, mfcc_delta2], axis=0)

# MFCC visualization
fig, axes = plt.subplots(3, 1, figsize=(12, 10))
librosa.display.specshow(mfcc_coefficients, sr=sr, hop_length=512,
                         x_axis='time', ax=axes[0])
axes[0].set_title('MFCC Coefficients')
axes[0].set_ylabel('MFCC')

librosa.display.specshow(mfcc_delta, sr=sr, hop_length=512,
                         x_axis='time', ax=axes[1])
axes[1].set_title('Delta MFCC')
axes[1].set_ylabel('Delta')

librosa.display.specshow(mfcc_delta2, sr=sr, hop_length=512,
                         x_axis='time', ax=axes[2])
axes[2].set_title('Delta-Delta MFCC')
axes[2].set_ylabel('Delta-Delta')

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

Noise Reduction Techniques

Noise reduction is a critical step in real-world applications. Different noise types require different techniques for optimal results:

Spectral subtraction estimates the ambient noise and subtracts it from the signal spectrum. It works well with stationary noise sources such as fans and air conditioning. A noise profile is extracted from silent segments where only noise is present, and this profile is subtracted from the entire signal.

Wiener filtering designs an optimal filter that maximizes the signal-to-noise ratio (SNR). It can be applied in both frequency and time domains, and adaptive versions can adjust to changing noise conditions.

def spectral_subtraction(noisy, noise_ref, fs, n_fft=2048):
    """Noise reduction via spectral subtraction."""
    # Estimate noise spectrum
    noise_fft = np.fft.rfft(noise_ref, n=n_fft)
    noise_power = np.mean(np.abs(noise_fft) ** 2)

    # Compute STFT of the signal
    f, t_stft, Zxx = scipy_signal.stft(noisy, fs=fs, nperseg=n_fft)

    # Subtract noise power (clamp negatives to zero)
    clean_power = np.maximum(np.abs(Zxx) ** 2 - noise_power, 0)
    clean_Zxx = np.sqrt(clean_power) * np.exp(1j * np.angle(Zxx))

    # Inverse STFT back to time domain
    _, clean_signal = scipy_signal.istft(clean_Zxx, fs=fs, nperseg=n_fft)
    return clean_signal

Underwater Acoustic Communication and Message Reconstruction

In my master's thesis, titled "ML-Based Reconstruction of Lost Acoustic Messages in Unmanned Underwater Vehicles," I worked on reconstructing acoustic messages that were lost or degraded during transmission through the underwater channel. This research addresses the unique challenges of underwater acoustic communication.

The underwater channel presents challenges fundamentally different from terrestrial communication:

• Multipath propagation: The signal reflects off the surface and seabed, arriving at the receiver at different times. These reflections overlap with the original signal, causing inter-symbol interference (ISI)
• Doppler shift: Relative movement between transmitter and receiver causes frequency shifts. Although underwater vehicle speeds are low, the speed of sound in water being much lower than in air makes the Doppler effect significant
• Ambient noise: Marine life, wave sounds, ship traffic, and thermal noise create a continuous noise floor
• Limited bandwidth: The usable frequency range underwater is extremely narrow, typically limited to a few kHz
• Distance-dependent attenuation: Acoustic signals in water suffer from both spreading and absorption losses as distance increases

In my thesis, I combined signal processing and machine learning techniques to reconstruct corrupted or completely lost acoustic messages under these challenging conditions. I built a comprehensive feature vector using FFT-based feature extraction for frequency-domain characteristics, spectrogram analysis for time-frequency domain features, and MFCC for perceptual features.

Deep learning models trained on these features achieved 90.8% command accuracy and 87.3% parameter recovery rate. Command accuracy refers to the rate at which transmitted acoustic commands were correctly identified; parameter recovery rate indicates how accurately the numerical parameters accompanying these commands (such as distance, angle, and speed) were reconstructed.

This work demonstrates the powerful results that can emerge when signal processing techniques are combined with artificial intelligence. Considering the uncertainties of the underwater environment and high noise levels, these success rates represent competitive results compared to existing work in the field.

Core Signal Processing Concepts

• Convolution: Combining two signals, the mathematical foundation of filtering
• Correlation: Measuring similarity between two signals
• Window Functions: Hamming, Hann, Blackman — to reduce spectral leakage
• Power Spectral Density (PSD): Distribution of signal power across frequencies
• Time-Frequency Analysis: Wavelet transform, STFT for time-varying frequency content

Application Domains

Signal processing techniques are widely used in:

• Audio and Music: Noise reduction, speech recognition, music information retrieval
• Telecommunications: Modulation, coding, channel equalization
• Medicine: ECG/EEG analysis, ultrasound imaging
• Radar and Sonar: Target detection, range measurement
• Seismology: Seismic data analysis, subsurface structure mapping

Summary

Signal processing and acoustic analysis play a foundational role across many engineering disciplines. FFT for frequency analysis, spectrograms for time-frequency visualization, and digital filters for noise removal are the essential tools in this field. The wavelet transform offers multi-resolution analysis beyond what the Fourier transform provides, while MFCC models the human auditory system to form the basis of speech recognition applications. The Python ecosystem — NumPy, SciPy, librosa, PyWavelets — makes it possible to apply these techniques quickly and effectively. Even in challenging domains like underwater acoustic communication, the combination of signal processing and AI techniques can achieve accuracy rates above 90%.