FUNCTION PERIODOGRAM,X,Y,D,S,WINDOW=window,WINFUN=w
;+
; PXY=PERIODOGRAM(X,Y,D,S)
; computes the cross power spectral density between
; sequences X and Y. If X and Y are the same sequence
; then a direct PSD is computed.
;
;
; INPUT PARAMETERS
; X		An array of data samples
; Y		An array of data samples
; D		The size of the analysis window
; S		The shift of the analysis window
;
; The window function is determined by the keyword WINDOW:
; WINDOW=0 (Default) All values equally weighted (no window)
; WINDOW=1 Hamming Window
; WINDOW=2 Bartlett Window
; WINDOW=3 Hann Window
; WINDOW=4 Truncated Gaussian
;
; The window function can be accessed via WINFUN keyword
;
; If X and Y are not the same length then the shorter
; one sets the length for the analysis.
;
; Example: Find the PSD of a sequence X at 500 points
; using a window of size 80 and shift of size 40.
;
; PSD=PERIODOGRAM(X,X,500,80,Window=1)
; PLOT,PSD
; Note that the PSD plot contains 500 points.
;
; H. Rhody  May, 1998
; Revised April 2000 to include window functions and
; remove extraneous parameter NF. Added access to
; window function via WINFUN.
;
; Reference:
;	S. L. Marple,
;	Digital Spectral Analysis with Applications,
;	Prentice Hall, 1987
; 	Appendix 5.C
;-

PSD=FLTARR(D)
XS=FLTARR(D)
YS=FLTARR(D)

IF NOT KEYWORD_SET(window) THEN window=0
k=(findgen(D)-(D-1)/2.)/(D-1)
CASE window OF
0:	w=FLTARR(D)+1 ; No window. Default
1:	BEGIN ; Hamming Window
		w=0.538+0.462*COS(2*!pi*k)
	END
2:	BEGIN  ; Bartlett Window
		w=1-2*ABS(k)
	END
3:	BEGIN ; Squared Cosine (Hann) Window
		w=0.5+0.5*cos(2*!pi*k)
	END
4:	BEGIN
		w=exp(-12.5*k^2)
	END
ELSE: MESSAGE, 'Unrecognized WINDOW selection'
ENDCASE
U=TOTAL(w)/D
w=w/U
NX=N_Elements(X)
NY=N_Elements(Y)
N=NX<NY
NP=(N-D)/S+1

FOR p=0,NP-1 DO BEGIN
	XS[0:D-1]=X(p*S:p*S+D-1)*w
	YS[0:D-1]=Y(p*S:p*S+D-1)*w
	XSF=FFT(XS)
	YSF=FFT(YS)
	PSD=PSD+XSF*CONJ(YSF)
ENDFOR
RETURN,ABS(PSD)/NP
END