; The baremusic procedure obtains the power spectrum from a string 
; of data. The string must be either a column or row vector.
; Use is as follows:
;                   baremusic,v,fs,avg,dtrnd,gxx,f,psisq
;    WHERE
;      |      v = variable name of data vector
; input|     fs = sampling frequency of data (hz)
;      |    avg = 1 to remove mean from data; 0 otherwise
;      |  dtrnd = 1 to detrend input data, 0 otherwise
;
;output|    gxx = one-sided spectrum
;      |      f = frequency vector corresponding to Gxx
;      |  psisq = mean sq value (if mean removed)
; No special window or smoothing used.
; This program developed from music.m written by B. Berger, 12/90. 
; Reference: Bendat & Piersol as well as D.E. Newland
; BWBerger 4/15/99
;***************************************************************************

pro baremusic,v,fs,avg,dtrnd,gxx,f,psisq

; Remove mean from time series if requested.

if dtrnd eq 1 then begin
   detren,v,v,dmn,tren,-999.
endif
if avg eq 1 then begin
   v=v-mean(v)
end

; Put data vector into column format.
rw=row(v)
cl=col(v)

if cl gt rw then begin
   v=transpose(v)
   rw=row(v)
   cl=col(v)
end
q=rw
; Find the number of points to use for requested fft type.
ev=2*floor(q/2)
if q gt ev then begin
   q=q-1
endif
n=q

; Do fft as requested.
a=fft(v(0:q-1),/double)
a=double(conj(a)*a)
t=(n-1)/double(fs)
s=2*t*a
gxx=s;
; Want only half of spectrum (remove frequency values larger than Nyquist).
half=n/2+1
gxx=transpose(gxx(0:half-2))
f=[dindgen(1,half-1)*fs/(n-1)]

psisq=(s(0)/(t*2))+(1/t)*total(s(1:half-2))

return
end