I'm creating a periodic time-domain signal by summing a series of sine waves with some amplitude, frequency, and phase. After performing a FFT with NumPy (reference), I construct a vector of the Fourier coefficients in the sine/cosine basis. Then I define a Fourier design matrix with sine/cosine elements evaluated at the different times and frequencies as shown below,
where T is the period of the signal, t are the time samples, and Nf is the number of Fourier modes used. I should be able to reconstruct the original signal by multiplying column vector of sine/cosine Fourier coefficients by the Fourier design matrix.
Here is the code I have to implement this:
import numpy as np
import matplotlib.pyplot as plt
# Function to create a real-valued signal with non-integer time samples
def create_signal(t, freq_components):
signal = np.zeros(len(t))
for amplitude, freq, phase in freq_components:
signal += amplitude * np.sin(2 * np.pi * freq * t + phase)
return signal
# Perform FFT and get the sine and cosine coefficients
def compute_fourier_coefficients(signal):
N = len(signal)
fft_result = np.fft.fft(signal)
# Initialize arrays for cosine and sine coefficients
cosine_coeffs = np.zeros(N // 2 + 1)
sine_coeffs = np.zeros(N // 2 + 1)
# Compute coefficients
for k in range(1, N // 2 + 1):
cosine_coeffs[k] = (2 / N) * fft_result[k].real
sine_coeffs[k] = -(2 / N) * fft_result[k].imag
# DC component (mean)
cosine_coeffs[0] = np.mean(signal)
return cosine_coeffs, sine_coeffs
# Create the Fourier design matrix with non-integer time samples
def create_fourier_design_matrix(t, num_modes, T):
N = len(t)
design_matrix = np.zeros((N, 2 * num_modes))
for k in range(1, num_modes + 1):
design_matrix[:, 2 * (k - 1)] = np.cos(2 * np.pi * k * t / T)
design_matrix[:, 2 * (k - 1) + 1] = np.sin(2 * np.pi * k * t / T)
return design_matrix
# Reconstruct the signal from the Fourier coefficients
def reconstruct_signal_from_design_matrix(fourier_design_matrix, cosine_coeffs, sine_coeffs):
num_modes = len(cosine_coeffs) - 1
coeffs = np.zeros(2 * num_modes)
coeffs[0::2] = cosine_coeffs[1:]
coeffs[1::2] = sine_coeffs[1:]
reconstructed_signal = fourier_design_matrix @ coeffs
reconstructed_signal += cosine_coeffs[0] # Add DC component to match original signal mean
return reconstructed_signal
# Parameters
N = 1024 # Number of samples
t = np.linspace(5000, 12000, N) # Non-integer time samples from 5000 to 12000
T = t[-1] - t[0] # Total duration
# Frequency components should correspond to actual frequencies in the signal
freq_components = [(1.0, 5 / T, 0), (0.5, 10 / T, np.pi/4), (0.2, 20 / T, np.pi/2)] # (amplitude, frequency, phase)
# Create the original signal
original_signal = create_signal(t, freq_components)
# Compute Fourier coefficients
cosine_coeffs, sine_coeffs = compute_fourier_coefficients(original_signal)
# Create Fourier design matrix
num_modes = N // 2
fourier_design_matrix = create_fourier_design_matrix(t, num_modes, T)
# Reconstruct the signal
reconstructed_signal = reconstruct_signal_from_design_matrix(fourier_design_matrix, cosine_coeffs, sine_coeffs)
# Plot the original and reconstructed signals
plt.plot(t, original_signal, label='Original Signal')
plt.plot(t, reconstructed_signal, label='Reconstructed Signal', linestyle='dashed')
plt.legend(loc='upper right')
plt.xlabel('time')
plt.ylabel('signal')
plt.show()
However, the reconstructed signal is offset in time from the original signal as shown below
What is going wrong?
You introduced the phase difference in the line:
t = np.linspace(5000, 12000, N)
Change this to:
t = np.linspace(0, 7000, N)
and you will be fine. The DFT is tacitly premised on the start coordinate being 0.