Как восстановить периодический сигнал из numpy.fft.fft, используя синусы, косинусы и матрицу проектирования Фурье?

Я создаю периодический сигнал во временной области, суммируя серию синусоидальных волн с некоторой амплитудой, частотой и фазой. После выполнения БПФ с помощью NumPy (ссылка) я строю вектор коэффициентов Фурье в базисе синуса/косинуса. Затем я определяю матрицу расчета Фурье с элементами синуса/косинуса, оцениваемыми в разное время и на разных частотах, как показано ниже:

где 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()

Однако восстановленный сигнал смещен во времени от исходного сигнала, как показано ниже.

Что происходит не так?

🤔 А знаете ли вы, что...
В Python можно легко работать с базами данных, такими как SQLite и MySQL.


1
68
1

Ответ:

Решено

Вы ввели разность фаз в строку:

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

Измените это на:

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

и с тобой все будет в порядке. ДПФ негласно предполагает, что начальная координата равна 0.