We consider the harmonic function of a field as follows:
x = amp* sin(2*pi*f*t)
Here, amp is the amplitude of the field, i.e, displacement or pressure, f is the frequency of the harmonic signal and the t is the time (parameter).
If we want to deal with RMS of the sinusoid signal than ampl/sqrt(2) should be used. Before performing Discrete Fourier Transform it is good practice to window signal in order to decreases leakage, etc.
Transforming the harmonic signal into frequency domain, two-sided spectrum is obtained. Normalization should be done done by dividing values of magnitude by energy of the window. If we to take signal "as it is" (which is equal to applying rectangular window) the spectrum is then divided by number of samples N. For different types of windows this factor that is equal to sum of all window samples.
Additionally we are only interested in one half of a spectrum, therefore amplitude of all samples must be multiplied by two to compensate the loss of energy, except of DC component which appears only once.
So the python code is as follows.
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal as mySignal
from scipy.signal import argrelextrema
plt.close('all')
f1 = 1
f2 = 2
f3 = 3
T1 = 1/f1
T2 = 1/f2
T3 = 1/f3
fs = f1*20
duration = 1/f3*8
npts = int(fs*duration)
t = np.arange(npts, dtype=float)/fs
amp1 = 1
amp2 = 2
amp3 = 3
signal = amp1 * np.sin((f1*duration) * np.linspace(0, 2*np.pi, npts)) + amp2 * np.sin((f2*duration) * np.linspace(0, 2*np.pi, npts)) + amp3 * np.sin((f3*duration) * np.linspace(0, 2*np.pi, npts))
signal3 =amp3 * np.sin((f3*duration) * np.linspace(0, 2*np.pi, npts))
plt.figure(1)
plt.subplot(2,1,1)
plt.plot(t, signal)
# Window signal
win_hamming = np.hamming(npts)
win_boxcar = mySignal.boxcar(npts)
win = win_boxcar
#win = win_hamming
signal = signal * win
plt.plot(t, signal)
plt.plot(t, signal3,linewidth=10)
plt.title("Boxcar window")
plt.ylabel("Amplitude")
plt.xlabel("Time,s")
sp = np.fft.fft(signal)
freq = np.fft.fftfreq(npts, 1.0/fs)
# Scale the magnitude of FFT by window energy and factor of 2,
# because we are using half of FFT.
sp_mag = np.abs(sp) * 2 / np.sum(win)
# To obtain RMS values, divide by sqrt(2)
sp_rms = sp_mag / np.sqrt(2)
# Shift both vectors to have DC at center
freq = np.fft.fftshift(freq)
sp_rms = np.fft.fftshift(sp_rms)
x = sp_rms*np.sqrt(2.0)
m = argrelextrema(x, np.greater) #array of indexes of the locals maxima
freq_local = [freq[m] for i in m]
sp_rms_local = [sp_rms[m]*np.sqrt(2.0) for i in m]
print(freq_local)
print(sp_rms_local)
plt.subplot(2,1,2)
plt.plot(freq, sp_rms*np.sqrt(2.0))
plt.plot(freq_local, sp_rms_local, 'rs')
plt.xlim( (0, f1*2) )
plt.grid('on')
plt.show()
