title: Can't find the right energy using scipy.signal.welch
body: <p>For a given discret time signal <code>x(t)</code> with spacing <code>dt</code> (which is equal to <code>1/fs</code>, <code>fs</code> being the sample rate), the energy is:</p>
  
  <pre><code>E[x(t)] = sum(abs(x)**2.0)/fs
  </code></pre>
  
  <p>Then I do a DFT of <code>x(t)</code>:</p>
  
  <pre><code>x_tf = np.fft.fftshift( np.fft.fft( x ) ) / ( fs * ( 2.0 * np.pi ) ** 0.5 )
  </code></pre>
  
  <p>and compute the energy again:</p>
  
  <pre><code>E[x_tf] = sum( abs( x_tf ) ** 2.0 ) * fs * 2 * np.pi / N
  </code></pre>
  
  <p>(here the factor <code>fs*2*np.pi/N</code> = pulsation spacing <code>dk</code>, the documentation of <code>fftfreq</code> gives more details about spacing in frequency domain), I have the same energy:</p>
  
  <pre><code>E[x(t)] = E[x_tf]
  </code></pre>
  
  <p>BUT... when I compute the power spectral density of <code>x(t)</code> using <code>scipy.signal.welch</code>, I can't find the right energy. <code>scipy.signal.welch</code> returns the vector of frequencies <code>f</code> and energy <code>Pxx</code> (or energy per frequency, depending on which <code>scaling</code> we enter in arguments of <code>scipy.signal.welch</code>).</p>
  
  <p>How can I find the same energy as <code>E[x(t)]</code> or <code>E[x_tf]</code> using <code>Pxx</code>? I tried to compute:</p>
  
  <pre><code>E_psd = sum(Pxx_den) / nperseg
  </code></pre>
  
  <p>where <code>nperseg</code> being the length of each segment of Welch algorithm, factors like <code>fs</code> and <code>np.sqrt(2*np.pi)</code> being cancelled out, and rescale E[x(t)] with <code>nperseg</code>, but without any success (orders of magnitude smaller than <code>E[x(t)]</code> )</p>
  
  <p>I used the following code to generate my signal:</p>
  
  <pre><code>#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)
  </code></pre>
  
  <p>and I did the following to get the power spectral density:</p>
  
  <pre><code>f, Pxx_den = signal.welch(x, fs )
  </code></pre>