# Compute the spectral density of a filtered version of the
# input spectrum.

spec.lfilter<-function(aneg,azero,apos,inputspec,inputfreq){
	nfreq<-length(inputfreq)
	nneg<-length(aneg)
	npos<-length(apos)
	
	neglags<-matrix(0,nfreq,nneg)
	for (k in 1:nneg) neglags[,k]<-aneg[k]*exp((1i)*2*pi*k*inputfreq)

	zerolag<-rep(azero,nfreq)
	
	poslags<-matrix(0,nfreq,npos)
	for (k in 1:npos) poslags[,k]<-apos[k]*exp((-1i)*2*pi*k*inputfreq)
		
	freqresponse<-apply(cbind(neglags,zerolag,poslags),1,sum)
	print(freqresponse)
	sqmodfreqresp<-Mod(freqresponse)^2
	outputspec<-sqmodfreqresp*inputspec
	speclfilter<-list(aneg=aneg,azero=azero,apos=apos,freq=inputfreq,spec=outputspec)
	plot(speclfilter$freq,speclfilter$spec,type="l")
	points(inputfreq,inputspec)
	return(speclfilter)
	}

inputfreq<-seq(0,.5,.01)
inputspec<-inputfreq*0+1

spec.lfilter(0,1,.5,inputspec,inputfreq)
spec.lfilter(c(0,0,.9),1,c(0,.9),inputspec,inputfreq)