1
$\begingroup$

I'm trying to check Parseval's theorm for Gaussian signal. It's well known that fourier transform of $\exp(-t^2)$ is $\sqrt{\pi}\exp(-\pi^2 k^2)$. So I implement it by using quad and simps. I think that the one is continuous integral, the other is integration which uses sample.

Code

from scipy import integrate
import numpy as np

def f(x):
    return (np.exp(-x**2))**2
def F(k):
    return ((np.pi**0.5)*np.exp((-np.pi**2)*(k**2)))**2

a=integrate.quad(f,-np.inf,np.inf)
print(a)
#>>>(1.2533141373155017, 4.4674503165883495e-09)
b=integrate.quad(F,-np.inf,np.inf)
print(b)
#>>>(1.2533141373155066, 1.592743224014847e-08)

a=b so above code follows Parseval's Theorm. The problem is on following code.

import scipy.fftpack as fft

N=50000
t  = np.linspace(-1000,1000, N)
h=np.exp(-t**2)
H=2*np.abs(fft.fftshift(fft.fft(h)/N))
freq=fft.fftshift(fft.fftfreq(H.shape[0],t[1]-t[0]))

S_h=integrate.simps(h**2,t)
print(S_h)
>>>1.25331413732

S_H=integrate.simps(H**2,freq)
print(S_H)
>>>1.25326400525e-06

S_h is not equals to S_H. What is the problem? How can I fix it??

$\endgroup$
  • $\begingroup$ try increasing N. and i might suggest making it a power of two which might make the FFT happier. $\endgroup$ – robert bristow-johnson Feb 14 at 4:02
  • $\begingroup$ If N>2500, S_H doesn't change regardless of increasing N.... $\endgroup$ – user10942748 Feb 14 at 6:30
  • $\begingroup$ your t still goes from -1000 to 1000. maybe that needs to go out as sqrt(N). your other other integration appears to be from -inf to inf (where does tjhat one stop?) $\endgroup$ – robert bristow-johnson Feb 14 at 7:41
  • $\begingroup$ and why is one integrate.quad() and the other is integrate.simps(). i know about simpson's rule (not particularly impressed with it). what does the quad do? $\endgroup$ – robert bristow-johnson Feb 14 at 7:48
  • $\begingroup$ try changing t = np.linspace(-1000,1000, N) to t = np.linspace(-(N**0.5),(N-1)**0.5, N) and let N be a power of 2 and get arbitrarily big. $\endgroup$ – robert bristow-johnson Feb 14 at 7:52
0
$\begingroup$

i dunno SciPy. but the two values appear to me very close. and when one does gaussians in the time or frequency domains, remember that once you sample, these single gaussians become replicated in the other domain and you have a pulse train of gaussians.

the Fourier Integral is not precisely the same as the approximation with finite limits and using Riemann summation, then you get the DFT (or fft).

With any discrete-time, discrete-frequency domain, the Riemann summation of an integral with finite limits might not be the same as the integral.

$\endgroup$
  • 1
    $\begingroup$ 1.253 is not "very close" to 1.253e-06 :) $\endgroup$ – endolith Feb 15 at 0:11
  • $\begingroup$ you're right, @endolith $\endgroup$ – robert bristow-johnson Feb 15 at 0:45
0
$\begingroup$

There is a simple idea behind Parseval's Theorem.
Unitary transform preserves the $ {L}_{2} $ Norm.
The Fourier Transform (All its variants) are indeed Unitary when the correct factor is in place.

The problem is with the way you sum the elements and factor them.
What you need to factor by $ \frac{1}{N} $ is the sum of squared absolute value of the series.

Have a look on the following MATLAB Code:

numSamples = 50000;
vT = linspace(-1000, 1000, numSamples);
vH = exp(-(vT .* vT));

vK = fft(vH);

% Summation by Parseval Theorem

sum(vH .* vH)
sum(vK .* conj(vK)) / numSamples

In MATLAB it yields

ans =

   31.3322

ans =

   31.3322

Why is that? Pay attention that the fft() applied by MATLAB and Python is not the Unitary form.
For the unitary form it should be factored by $ \frac{1}{\sqrt{N}} $.

Yet in Parseval we see:

$$ \sum_{n = 0}^{N - 1} {x}^{2} \left[ n \right] = \frac{1}{N} \sum_{k = 0}^{N - 1} {\left| X \left[ k \right] \right|}^{2} $$

The factor $ \frac{1}{N} $ is exactly $ {\left( \frac{1}{\sqrt{N}} \right)}^{2}$.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.