Теорема Парсеваля с Numpy FFT не выполняется

avatar
Alberto
3 мая 2021 в 12:51
214
1
2

Я пытаюсь определить общую энергию, зарегистрированную детектором во временной области, с помощью его спектра. Первым шагом после выполнения быстрого преобразования Фурье с библиотекой БПФ Numpy было подтверждение теоремы Парсеваля.

Согласно теореме, полная энергия во временной области и в частотной области должна быть одинаковой. У меня есть две проблемы, которые я не могу решить.

  • Я могу подтвердить теорему, если не использую правильные единицы измерения для оси X во время интеграции np.trapz(). Как только я использую свои фактические точки/частоты выборки, результат отключается. Я не понимаю, почему это так, и мне интересно, могу ли я применить нормализацию для устранения этой ошибки.
  • Я не могу подтвердить теорему, когда применяю смещение постоянного тока к сигналу (раскомментируйте строку f = np.sin(np.pi**t)*).

Ниже мой код с примером функции синуса.

# Python code

import matplotlib.pyplot as plt
import numpy as np

# Create a Sine function
dt = 0.001 # Time steps
t = np.arange(0,10,dt) # Time array
f = np.sin(np.pi*t) # Sine function
# f = np.sin(np.pi*t)+1 # Sine function with DC offset
N = len(t) # Number of samples

# Energy of function in time domain
energy_t = np.trapz(abs(f)**2)

# Energy of function in frequency domain
FFT = np.sqrt(2) * np.fft.rfft(f) # only positive frequencies; correct magnitude due to discarding of negative frequencies
FFT[0] /= np.sqrt(2) # DC magnitude does not have to be corrected
FFT[-1] /= np.sqrt(2) # Nyquist frequency does not have to be corrected
frq = np.fft.rfftfreq(N,d=dt) # FFT frequenices

# Energy of function in frequency domain
energy_f = np.trapz(abs(FFT)**2) / N

print('Parsevals theorem fulfilled: ' + str(energy_t - energy_f))

# Parsevals theorem with proper sample points

energy_t = np.trapz(abs(f)**2, x=t)
energy_f = np.trapz(abs(FFT)**2, x=frq) / N

print('Parsevals theorem NOT fulfilled: ' + str(energy_t - energy_f))
Источник
Matt Timmermans
3 мая 2021 в 13:08
0

trapz придает только половину веса первому и последнему образцам. Используйте sum или добавьте первый образец в конец, так как ожидается, что и частоты, и сигнал будут циклическими.

Alberto
4 мая 2021 в 08:55
0

Оба метода сработали, спасибо!

Ответы (1)

avatar
Cris Luengo
3 мая 2021 в 14:53
3

БПФ вычисляет дискретное преобразование Фурье (ДПФ), которое не совпадает с преобразованием Фурье (в непрерывной области).

Для ДПФ теорема Парсеваля утверждает, что сумма квадратов амплитуд дискретного сигнала равна сумме квадратов амплитуд ДПФ сигнала. Интеграция не требуется, поэтому не следует использовать trapz. Просто используйте sum.


Note that a discrete signal is a set of samples x[n] at n=0..N-1 . Анализ Фурье в дискретной области и все связанные операции учитывают только n, а не t. Частота выборки и фактическое время записи этих выборок не имеет значения для этих анализов. Аналогично, ДПФ создает набор отсчетов X[k] в k=0..=0..=0..=0..<14064049651587>n конкретные f или ω, относящиеся к любой частоте дискретизации.

Теперь можно связать n с t, потому что мы знаем частоту дискретизации, и можно связать k<<140640496560>60490<140640496560> > потому что мы знаем частоту дискретизации. Но эти преобразования не должны заставлять нас думать, что X[k] является выборкой преобразования Фурье в непрерывной области исходного сигнала в непрерывной области. И они особо не должны наводить нас на мысль, что мы можем интерполировать X[k].

Reconstructing the samples x[n] is accomplished by adding N sinusoids with parameters given by X[k ]. «Между» этими компонентами DFT не должно быть ничего. Их интерполяция означала бы добавление синусоид, которых нет в выборках x[n].

trapz использует линейную интерполяцию для получения оценки интеграла и поэтому не подходит для использования в дискретном анализе Фурье.

Alberto
4 мая 2021 в 08:55
0

Да, sum выполняет свою работу, и я получаю правильные результаты, спасибо! Однако, поскольку полная энергия сохраняется для fft, почему интегралы во временной области по времени и в частотной области по частоте не равны? @Крис

Cris Luengo
4 мая 2021 в 13:25
0

@Alberto: использование trapz подразумевает линейную интерполяцию. Я добавил обсуждение к ответу, почему это неуместно.