Hello matplotlib users.
I'm new to signal processing and I've read that RMS could be found from
a PSD. I'm interested in as I would further like to know energy in a
signal through it's frequencies.
My problem is I don't find how to calculate the RMS from the PSD output.
It seems it's a matter of scale (frequencies bandwith is taken in
account already).
I wrote a test case with a simple sinus. I should be able to find the
same RMS value from the PSD method and direct RMS over signal method.
Could you please have a look and tell me how to find good RMS value from
PSD output?
Thanks
#!/usr/bin/python
# * coding: utf8 *
import matplotlib, platform
if platform.system() == 'Linux' :
matplotlib.use("gtk")
import pylab
import scipy
## PSD vs RMS
#Parameters
samplerate = 48000
nfft = 1024*2
graph = False
#create 1 sec sinus signal
t = scipy.arange(0, 1 , 1/float(samplerate))
signal = .25*scipy.sin(2*scipy.pi*(samplerate/10.)*t)
print len(signal)
#RMS of an array
def RMS(data):
rms = data**2
rms = scipy.sqrt(rms.sum()/len(data))
return rms
#PSD of an array. I want this to return the RMS
def RMSfromPSD(data) :
y, x = pylab.psd(data, NFFT = nfft, Fs = samplerate)
##Calculate the RMS
#The energy returned by PSD depends on FFT size
freqbandwith = x[1]
y = y*freqbandwith
#The energy returned by PSD depends on Samplerate
y = y/float(samplerate)
#Summing the power in freq domain to get RMS
rms = scipy.sqrt(y.sum())
return rms
print "RMS method", RMS(signal)
print "RMS using PSD method", RMSfromPSD(signal)
#Graph
if graph == True :
pylab.subplot(211)
pylab.plot(t,signal)
pylab.subplot(212)
pylab.psd(signal, nfft, samplerate)
pylab.show()
