Teorema di Parseval in Python

Question

Sto cercando di ottenere un po 'di controllo sulla funzionalità fft di Python, e una delle cose strane su cui sono incappato è che non sembra applicarsi lo Parseval's theorem, dato che dà una differenza di circa 50 ora, mentre dovrebbe essere 0.

import numpy as np 
import matplotlib.pyplot as plt 
import scipy.fftpack as fftpack 

pi = np.pi 

tdata = np.arange(5999.)/300 
dt = tdata[1]-tdata[0] 

datay = np.sin(pi*tdata)+2*np.sin(pi*2*tdata) 
N = len(datay) 

fouriery = abs(fftpack.rfft(datay))/N 

freqs = fftpack.rfftfreq(len(datay), d=(tdata[1]-tdata[0])) 

df = freqs[1] - freqs[0] 

parceval = sum(datay**2)*dt - sum(fouriery**2)*df 
print parceval 

plt.plot(freqs, fouriery, 'b-') 
plt.xlim(0,3) 
plt.show() 

sono abbastanza sicuro che si tratta di un fattore di normalizzazione, ma non mi sembra di essere in grado di trovare, come di tutte le informazioni che posso trovare su questa funzione è il scipy.fftpack.rfft documentation.

Answer 1

Il fattore di normalizzazione proviene dal tentativo di applicare il teorema di Parseval per la trasformata di Fourier di un segnale continuo a una sequenza discreta. Sul pannello laterale di the wikipedia article on the Discrete Fourier transform c'è qualche discussione sulla relazione tra la trasformata di Fourier, la serie di Fourier, la Trasformata di Fourier discreta e il campionamento con Dirac combs.

per fare una lunga storia breve, Parseval's theorem, when applied to DFTs, non richiede l'integrazione, ma sommatoria: un 2*pi si sta creando da multipliying da dt e df tuoi sommatorie.

Si noti inoltre che, poiché si utilizza scipy.fftpack.rfft, ciò che si ottiene non è propriamente la DFT dei dati, ma solo la metà positiva di esso, poiché il negativo sarebbe simmetrico ad esso. Quindi, dato che stai aggiungendo solo metà dei dati, più il 0 nel termine DC, manca lo 2 per arrivare allo 4*pi trovato da @unutbu.

In ogni caso, se datay tiene la sequenza, è possibile verificare il teorema di Parseval come segue:

fouriery = fftpack.rfft(datay) 
N = len(datay) 
parseval_1 = np.sum(datay**2) 
parseval_2 = (fouriery[0]**2 + 2 * np.sum(fouriery[1:]**2))/N 
print parseval_1 - parseval_2 

Utilizzando scipy.fftpack.fft o numpy.fft.fft la seconda sommatoria non ha bisogno di assumere una forma così strana:

fouriery_1 = fftpack.fft(datay) 
fouriery_2 = np.fft.fft(datay) 
N = len(datay) 
parseval_1 = np.sum(datay**2) 
parseval_2_1 = np.sum(np.abs(fouriery_1)**2)/N 
parseval_2_2 = np.sum(np.abs(fouriery_2)**2)/N 
print parseval_1 - parseval_2_1 
print parseval_1 - parseval_2_2