How to reconstruct periodic signal from numpy.fft.fft using sines, cosines, and a Fourier design matrix?

Question

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?

Answer 1

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.