0
$\begingroup$

I have the following code which generates a single simple sine wave with amplitude $1$ and then calculates the Fourier transform, normalises it, and finally takes the sum of all the components.

However I've observed something strange, if I generate a $10\textrm{ Hz}$ wave and a $100\textrm{ Hz}$ wave. The amplitudes are pretty damn similar as expected and they are 0.9999483519908307 and 0.9997854546246567 respectively.

However the total sum of all the components after normalisation are 1.01723486798 and 1.17263231047 respectively, which is approximately a $15\%$ difference! Why is this occurring?

Below are the plots of the spectrums and the python code below.

10Hz

100Hz

mport matplotlib.pyplot as plt
import numpy as np
from decimal import Decimal

#number of samples
iNumSamples=10000

#period of samples (i.e. sampling at 10Khz)
Ts=1/iNumSamples

#time axis used to sample 10,000 times for a second
t = np.linspace(0, 1, iNumSamples)

#10Hz Signal
signal1=np.sin(t*20*np.pi)

#100Hz Signal
signal2=np.sin(t*200*np.pi)


#creates a new figure object 
freqFig=plt.figure(2)
#change title
freqFig.canvas.set_window_title('Frequency domain') 

#fast fourier transform of 10Hz signal
absoluteValues1=abs(np.fft.rfft(signal1)) 

#fast fourier transform of 100Hz
absoluteValues2=abs(np.fft.rfft(signal2)) 

#find the frequencies of the rtansform
frequencies=np.fft.rfftfreq(iNumSamples, d=Ts)


##Normalises (rescales) the fft in a sensible way so that the corresponding freq component of a sine wave is 1 if the sine wave amplitude is 1
##non zero freq amplitudes (even when using rfft) are halfed as there are negative and positive components, hence first multiple by 2

#first create empty array of correct size
normalisedValues1=np.zeros(absoluteValues1.size, dtype=Decimal)
normalisedValues2=np.zeros(absoluteValues2.size, dtype=Decimal)

#normalise 0 component
normalisedValues1[0]=absoluteValues1[0]/iNumSamples 
normalisedValues2[0]=absoluteValues2[0]/iNumSamples 
#normalise non zero components
normalisedValues1[1:]=absoluteValues1[1:]*(2/iNumSamples)
normalisedValues2[1:]=absoluteValues2[1:]*(2/iNumSamples)

print("peak value signal 1:")
print(np.amax(normalisedValues1))
print("peak value signal 2:")
print(np.amax(normalisedValues2))

print("total sum of all componenets in spectrum of 10 Hz")
print(np.sum(normalisedValues1))

print("total sum of all componenets in spectrum of 100 Hz")
print(np.sum(normalisedValues2))

#creates a new figure object 
freqFig=plt.figure(1)
#change title
freqFig.canvas.set_window_title('Frequency domain') 
plt.plot(frequencies,normalisedValues1)

#creates a new figure object 
freqFig=plt.figure(2)
#change title
freqFig.canvas.set_window_title('Frequency domain') 
plt.plot(frequencies,normalisedValues2)

plt.show()
$\endgroup$
0
$\begingroup$

Below is the corrected source code which outputs values which are more as you would expect. I've marked the parts of the source I have corrected.

First, you did not draw a full period of the signal. You had

t = np.linspace(0, 1, iNumSamples)

which I have corrected to

t = np.linspace(0, 1, iNumSamples, endpoint=False)

For example, with iNumSamples=10 your code would generate array([ 0,0.11111111,0.22222222,0.33333333,0.44444444,0.55555556,0.66666667,0.77777778,0.88888889,1.]). But, what you want to have is t=[0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9], which is what my line is doing. If the time would go to up to 1, you would already start with the first sample of the following period.

The second point is a wrong understanding of the Parsevals Theorem. The theorem states

$$\|x[n]\|^2=\|X[k]\|^2$$

where $x[n]$ is the time domain sequence and $X[k]$ is the corresponding Fourier transform (Here, I've assumed unitary Fourier transform, otherwise some normalization factor would occur). The main thing here is that

$$\|X[k]\|^2=\sum_k|X[k]|^2$$

and not

$$\|X[k]\|^2\neq \sum_k|X[k]|$$

which is what you have calculated. So, you need to sum up the squares of the absolute values to calculate the Energy in the frequency domain, and this will match the energy in the time domain.

Finally, here's the source:

import matplotlib.pyplot as plt
import numpy as np
from decimal import Decimal

#number of samples
iNumSamples=10000

#period of samples (i.e. sampling at 10Khz)
Ts=1/iNumSamples

#time axis used to sample 10,000 times for a second
t = np.linspace(0, 1, iNumSamples, endpoint=False)  ##### Corrected

#10Hz Signal
signal1=np.sin(t*20*np.pi)

#100Hz Signal
signal2=np.sin(t*200*np.pi)


#creates a new figure object 
freqFig=plt.figure(2)
#change title
freqFig.canvas.set_window_title('Frequency domain') 

#fast fourier transform of 10Hz signal
absoluteValues1=abs(np.fft.rfft(signal1)) 

#fast fourier transform of 100Hz
absoluteValues2=abs(np.fft.rfft(signal2)) 

#find the frequencies of the rtansform
frequencies=np.fft.rfftfreq(iNumSamples, d=Ts)


##Normalises (rescales) the fft in a sensible way so that the corresponding freq component of a sine wave is 1 if the sine wave amplitude is 1
##non zero freq amplitudes (even when using rfft) are halfed as there are negative and positive components, hence first multiple by 2

#first create empty array of correct size
normalisedValues1=np.zeros(absoluteValues1.size, dtype=Decimal)
normalisedValues2=np.zeros(absoluteValues2.size, dtype=Decimal)

#normalise 0 component
normalisedValues1[0]=absoluteValues1[0]/iNumSamples 
normalisedValues2[0]=absoluteValues2[0]/iNumSamples 
#normalise non zero components
normalisedValues1[1:]=absoluteValues1[1:]*(2/iNumSamples)
normalisedValues2[1:]=absoluteValues2[1:]*(2/iNumSamples)

print("peak value signal 1:")
print(np.amax(normalisedValues1))
print("peak value signal 2:")
print(np.amax(normalisedValues2))

print("total sum of all componenets in spectrum of 10 Hz")
print(np.sum(normalisedValues1**2))  #### Corrected

print("total sum of all componenets in spectrum of 100 Hz")
print(np.sum(normalisedValues2**2))   #### Corrected

#creates a new figure object 
freqFig=plt.figure(1)
#change title
freqFig.canvas.set_window_title('Frequency domain') 
plt.plot(frequencies,normalisedValues1)
plt.xlim((0, 300))

#creates a new figure object 
freqFig=plt.figure(2)
#change title
freqFig.canvas.set_window_title('Frequency domain') 
plt.plot(frequencies,normalisedValues2)
plt.xlim((0,300))

Program output:

peak value signal 1:
1.0
peak value signal 2:
1.0
total sum of all componenets in spectrum of 10 Hz
1.0
total sum of all componenets in spectrum of 100 Hz
1.0
$\endgroup$
1
$\begingroup$

That's a subtle sampling issues and you will get your desired result by replacing this line:

t = np.linspace(0, 1, iNumSamples)

with the following:

t = np.linspace(0, 1-Ts, iNumSamples)

That being said, it's really a virtual sucess and you would never rely on such perfectly matching numbers in practice...

$\endgroup$
0
$\begingroup$

The sum of the squared amplitudes will give you the energy of the signal (Parseval). And this should be the same for both signals.

The sum of the (unsquared) amplitudes does not make much sense.

$\endgroup$
0
$\begingroup$

It doesn't. The fact that you changed the scale factor between the different frequency plots is what changed the sum. Try using the same normalization constant for both plots.

The reason your scale factors changed is that the magnitude of the frequency peak of a sinusoid is not the same as the magnitude of the nearest FFT result bin. Look up scalloping loss. So your plots were normalized by two different slightly bogus scale factors.

$\endgroup$
  • $\begingroup$ I did use the same normalisation constant for both plots. Its the one described in this link which provides a frequency component of amplitude 1 when the corresponding sine wave has a frequency of 1. matteolandi.blogspot.co.uk/2010/08/… $\endgroup$ – generic purple turtle Mar 6 '17 at 11:33

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.