Parseval's Theorem does not hold for FFT of a

Question

Thanks in advance for any help on this subject. I've recently been trying to work out Parseval's theorem for discrete fourier transforms when noise is included. I based my code from this code.

What I expected to see is that (as when no noise is included) the total power in the frequency domain is half that of the total power in the time-domain, as I have cut off the negative frequencies.

However, as more noise is added to the time-domain signal, the total power of the fourier transform of the signal+noise becomes much less than half of the total power of the signal+noise.

My code is as follows:

import numpy as np
import numpy.fft as nf
import matplotlib.pyplot as plt

def findingdifference(randomvalues):  

    n               = int(1e7) #number of points
    tmax            = 40e-3 #measurement time
    f1              = 30e6 #beat frequency

    t               = np.linspace(-tmax,tmax,num=n) #define time axis
    dt              = t[1]-t[0] #time spacing

    gt              = np.sin(2*np.pi*f1*t)+randomvalues #make a sin + noise

    fftfreq         = nf.fftfreq(n,dt) #defining frequency (x) axis
    hkk             = nf.fft(gt) # fourier transform of sinusoid + noise
    hkn             = nf.fft(randomvalues) #fourier transform of just noise

    fftfreq         = fftfreq[fftfreq>0] #only taking positive frequencies
    hkk             = hkk[fftfreq>0]
    hkn             = hkn[fftfreq>0]

    timedomain_p    = sum(abs(gt)**2.0)*dt #parseval's theorem for time
    freqdomain_p    = sum(abs(hkk)**2.0)*dt/n # parseval's therom for frequency
    difference      = (timedomain_p-freqdomain_p)/timedomain_p*100 #percentage diff

    tdomain_pn   = sum(abs(randomvalues)**2.0)*dt #parseval's for time
    fdomain_pn   = sum(abs(hkn)**2.0)*dt/n # parseval's for frequency
    difference_n    = (tdomain_pn-fdomain_pn)/tdomain_pn*100 #percent diff

    return difference,difference_n

def definingvalues(max_amp,length):

    noise_amplitude     = np.linspace(0,max_amp,length) #defining noise amplitude
    difference          = np.zeros((2,len(noise_amplitude)))
    randomvals          = np.random.random(int(1e7)) #defining noise

    for i in range(len(noise_amplitude)):
        difference[:,i] = (findingdifference(noise_amplitude[i]*randomvals))

    return noise_amplitude,difference

def figure(max_amp,length):

    noise_amplitude,difference = definingvalues(max_amp,length)

    plt.figure()
    plt.plot(noise_amplitude,difference[0,:],color='red')
    plt.plot(noise_amplitude,difference[1,:],color='blue')
    plt.xlabel('Noise_Variable')
    plt.ylabel(r'Difference in $\%$')
    plt.show()

    return
figure(max_amp=3,length=21)

My final graph looks like this figure. Am I doing something wrong when working this out? Is there an physical reason that this trend occurs with added noise? Is it to do with doing a fourier transform on a not perfectly sinusoidal signal? The reason I am doing this is to understand a very noisy sinusoidal signal that I have real data for.

Answer 1

Parseval's theorem holds in general if you use the whole spectrum (positive and negative) frequencies to compute the power.

The reason for the discrepancy is the DC (f=0) component, which is treated somewhat special.

First, where does the DC component come from? You use np.random.random to generate random values between 0 and 1. So on average you raise the signal by 0.5*noise_amplitude, which entails a lot of power. This power is correctly computed in the time domain.

However, in the frequency domain, there is only a single FFT bin that corresponds to f=0. The power of all other frequencies is distributed over two bins, only the DC power is contained in a single bin.

By scaling the noise you add DC power. By removing the negative frequencies you remove half the signal power, but most of the noise power is located in the DC component which is used fully.

You have several options:

1. Use all frequencies to compute the power.
2. Use noise without a DC component: randomvals = np.random.random(int(1e7)) - 0.5
3. "Fix" the power calculation by removing half of the DC power: hkk[fftfreq==0] /= np.sqrt(2)

I'd go with option 1. The second might be OK and I don't really recommend 3.

Finally, there is a minor problem with the code:

fftfreq         = fftfreq[fftfreq>0] #only taking positive frequencies
hkk             = hkk[fftfreq>0]
hkn             = hkn[fftfreq>0]

This does not really make sense. Better change it to

hkk             = hkk[fftfreq>=0]
hkn             = hkn[fftfreq>=0]

or completely remove it for option 1.