R语言 功率谱密度分析

library(psd)
fs=1000;
nfft = 1024;   
t=seq(0,1,1/fs)
xn=3*cos(2*pi*100*t)+3*cos(2*pi*300*t)+10*cos(2*pi*200*t)
#par(mfrow=c(1,2))
#plot(xn,t='l')
psdc=pspectrum(xn,x.frqsamp=fs,Nyquist.normalize=FALSE,plot=FALSE)
#c(psdc)
plot(psdc,log="dB", main="power spectral density estimates",lwd=2, ci.col=NA ,ylim=c(0,50),yaxs="i",xlim=c(0,fs/2))
sumpsd=sum(psdc$spec)
percent=function(x1,x2)
{sum(psdc$spec[x1:x2])/sumpsd}
a=percent(180,220)
c(a)#0.84723