如何使用Python从已知时延的两个时移信号叠加结果中重构原始信号S(t)

Great question! Reconstructing the original signal ( S(t) ) from the delayed sum ( S_{out}(t) = S(t) + S(t-dt) ) is a classic linear signal processing problem, and there are two solid approaches to implement this in Python. Let's walk through both with code examples.

Method 1: Frequency-Domain Inverse Filtering (Offline Processing) #

This is the most direct method, leveraging Fourier transforms to turn the time-domain superposition into a simple division in the frequency domain.

How it works #

Take the Fourier transform of both sides of the equation:

( \mathcal{F}{S_{out}(t)} = \mathcal{F}{S(t)} + \mathcal{F}{S(t-dt)} )

Using the time-shifting property of Fourier transforms (( \mathcal{F}{S(t-dt)} = \mathcal{F}{S(t)}e^{-j\omega dt} )), we get:

( F_{out}(\omega) = F_S(\omega) \left(1 + e^{-j\omega dt}\right) )

Rearranged to solve for the original signal's Fourier transform:

( F_S(\omega) = \frac{F_{out}(\omega)}{1 + e^{-j\omega dt}} )

Then we just take the inverse Fourier transform of ( F_S(\omega) ) to get back ( S(t) ).

Critical Note: At frequencies where ( 1 + e^{-j\omega dt} = 0 ) (i.e., ( \omega dt = \pi + 2k\pi ), or ( f = (2k+1)/(2dt) )), the denominator is zero, which would cause numerical explosions. We fix this by adding a tiny epsilon value to the denominator to regularize it.

Python Code Example #

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

# Step 1: Generate sample data
fs = 1000  # Sampling frequency (Hz)
dt = 0.01  # Known time delay (seconds)
t = np.linspace(0, 1, fs, endpoint=False)  # Time vector

# Original signal (example: a mix of sine waves)
S = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*20*t)

# Create the delayed sum signal
S_out = S + np.roll(S, int(dt*fs))  # np.roll shifts the signal by dt*fs samples

# Step 2: Frequency-domain reconstruction
N = len(S_out)
freqs = fftfreq(N, 1/fs)

# Compute Fourier transforms
F_out = fft(S_out)
omega = 2 * np.pi * freqs

# Compute the filter kernel, with regularization to avoid division by zero
epsilon = 1e-6
filter_kernel = 1 / (1 + np.exp(-1j * omega * dt) + epsilon)

# Apply filter and inverse transform
F_S = F_out * filter_kernel
S_recon = np.real(ifft(F_S))  # Take real part to discard tiny imaginary artifacts

# Step 3: Plot results
plt.figure(figsize=(12, 6))
plt.subplot(2,1,1)
plt.plot(t, S, label='Original S(t)')
plt.plot(t, S_out, label='S_out(t) = S(t) + S(t-dt)', alpha=0.7)
plt.title('Original and Mixed Signal')
plt.legend()

plt.subplot(2,1,2)
plt.plot(t, S, label='Original S(t)')
plt.plot(t, S_recon, label='Reconstructed S(t)', linestyle='--')
plt.title('Reconstruction via Frequency-Domain Filtering')
plt.legend()
plt.tight_layout()
plt.show()

Method 2: Time-Domain Iterative Reconstruction (Real-Time Friendly) #

If you need a method that doesn't rely on Fourier transforms (e.g., for real-time or streaming processing), iterative subtraction works well. The idea is to start with an initial guess of ( S(t) ), then repeatedly subtract the delayed version of your guess from ( S_{out}(t) ) to refine the estimate.

How it works #

1. Start with an initial guess: ( S_0(t) = S_{out}(t) )
2. Iterate using: ( S_{n+1}(t) = S_{out}(t) - S_n(t-dt) )
3. After a few iterations, ( S_n(t) ) will converge to the original ( S(t) ) (since each iteration removes the delayed component from the previous guess)

Python Code Example #

import numpy as np
import matplotlib.pyplot as plt

# Step 1: Generate sample data (same as before)
fs = 1000
dt = 0.01
t = np.linspace(0, 1, fs, endpoint=False)
S = np.sin(2*np.pi*5*t) + 0.5*np.sin(2*np.pi*20*t)
S_out = S + np.roll(S, int(dt*fs))

# Step 2: Iterative reconstruction
num_iterations = 10
S_recon = S_out.copy()  # Initial guess

for _ in range(num_iterations):
    # Subtract the delayed version of the current reconstruction
    S_recon = S_out - np.roll(S_recon, int(dt*fs))

# Step 3: Plot results
plt.figure(figsize=(12, 6))
plt.subplot(2,1,1)
plt.plot(t, S, label='Original S(t)')
plt.plot(t, S_out, label='S_out(t) = S(t) + S(t-dt)', alpha=0.7)
plt.title('Original and Mixed Signal')
plt.legend()

plt.subplot(2,1,2)
plt.plot(t, S, label='Original S(t)')
plt.plot(t, S_recon, label='Reconstructed S(t) (10 iterations)', linestyle='--')
plt.title('Reconstruction via Time-Domain Iteration')
plt.legend()
plt.tight_layout()
plt.show()

Which Method to Choose? #

• Frequency-Domain Filtering: Fast, efficient for offline processing. Best when you have the entire signal available upfront. Just remember to handle the division-by-zero singularities with regularization.
• Time-Domain Iteration: No Fourier transforms needed, making it suitable for real-time or streaming scenarios. Converges quickly (usually 5-10 iterations are enough for most cases).

内容的提问来源于stack exchange，提问作者Mariusz