Для данного дискретного сигнала времени x(t)
с интервалом dt
(который равен 1/fs
, fs
- частота дискретизации), энергия равна:
E[x(t)] = sum(abs(x)**2.0)/fs
Тогда I сделать ДПФ из x(t)
:
x_tf = np.fft.fftshift( np.fft.fft( x ) ) / ( fs * ( 2.0 * np.pi ) ** 0.5 )
и снова вычислить энергию:
E[x_tf] = sum( abs( x_tf ) ** 2.0 ) * fs * 2 * np.pi / N
(здесь коэффициент fs*2*np.pi/N
= расстояние между пульсациями dk
, документация fftfreq
дает более подробную информацию о разносе в частотной области), у меня та же энергия:
E[x(t)] = E[x_tf]
НО ... когда я вычисляю спектральную плотность мощности x(t)
с использованием scipy.signal.welch
, я могу ' найти правильную энергию. scipy.signal.welch
возвращает вектор частот f
и энергию Pxx
(или энергию на частоту, в зависимости от того, какое scaling
мы введем в аргументах scipy.signal.welch
).
Как я могу найти ту же энергию, что и E[x(t)]
или E[x_tf]
, используя Pxx
? Я попытался вычислить:
E_psd = sum(Pxx_den) / nperseg
где nperseg
- длина каждого сегмента алгоритма Уэлча, факторы как fs
и np.sqrt(2*np.pi)
отменены, и масштабирование E [x (t)] с nperseg
, но безуспешно (на порядки меньше, чем E[x(t)]
)
Я использовал следующий код для генерации моего сигнала:
#Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by 0.001 V**2/Hz of white noise sampled at 10 kHz.
fs = 10e3 #sampling rate, dt = 1/fs
N = 1e5
amp = 2*np.sqrt(2)
freq = 1234.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp*np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
и сделал следующее для получить спектральную плотность мощности:
f, Pxx_den = signal.welch(x, fs )