如何使用正弦、余弦和傅里叶设计矩阵从 numpy.fft.fft 重建周期信号？

Question

我通过将一系列具有一定幅度、频率和相位的正弦波相加来创建周期性时域信号。使用 NumPy（reference）执行 FFT 后，我在正弦/余弦基础上构建了傅里叶系数的向量。然后，我定义一个傅立叶设计矩阵，其中包含在不同时间和频率评估的正弦/余弦元素，如下所示，

其中 T 是信号的周期，t 是时间样本，Nf 是使用的傅里叶模式的数量。我应该能够通过将正弦/余弦傅里叶系数的列向量乘以傅里叶设计矩阵来重建原始信号。

这是我必须实现的代码：

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()

但是，重建信号与原始信号存在时间偏移，如下所示

出了什么问题？

Answer 1

您自己介绍了相位差！

换线

t = np.linspace(5000, 12000, N)

到

t = np.linspace(0, 7000, N)

你会没事的。 DFT默认起始坐标为0。