function [s100,F,y] = multispec(x,sr,ndt,txxu,txv,txu,txx);
%
%function [s100,F,y] = multispec(x,sr,ndt,txxu,txv,txu,txx);
%
% 13 Jan 2004
% David S. Ciochetto
%
% A routine to interface to Matlab's psd.m power spectrum estimate. It correctly normalizes the output
% from psd such that the integral of the area under the power spectra equals the variance of the original
% signal. It is set up to use the power spectral estimate from a number of ensembles of a signal (x)
% and to take the average to reduce the variance of the noise in frequency space (plotted against frequency).
%
% The method does not use overlapping internally so an estimate of the 95% confidence bounds can be added.
%
% It will automatically make a plot of the estimate of the power spectrum. Optionally, you can get the data
% for the (x,y) data to reconstruct the plot. It can also return the upper and lower 95% confidence bounds.
%
% Input: 
%  x   = matrix of ensembles of data in columns, one ensemble per column.
%  sr  = the sample rate, default units are Hz but can be nay units specified in txxu for proper labeling.
%  ndt = if you want to remove linear trends from the data ('y' or 'n').
%     Text to label the plot it generates with. (all text strings)
%  txxu = the units of the frequency axis (e.g. 'Hz' = 1 sample / second, could use 1 sample / day)
%  txv  = name of variable (used for ylabel), default is 'x(t)'
%  txu  = units of variable (used for ylabel), default is 'units of x'
%  txx  = extra text to add after the title to identify the data, default is 'sample data'
%
% Output:
%  s100 = estimate of the power spectrum. (the result or the dependent variable for a plot PSD(f) v. f.
%  F    = frequencies at which the estimates were made (frequency axis data or independent variable)
%  y    = 95% confidence limits. Matrix with two columns, y(:,1) is upper bound for 95% confidence, y(2) is the lower.
%
%
% If nothing is specified for ndt, it removes a linear trend from each ensemble. It is highly suggested to remove
% trends from the data to get rid of low frequency masking of higher frequency signals of interest.
%
% Based on my spectral analysis file ts001.m that I developed for Alex Hays time series class. 
%
% It will return the spectrum averaged over all of the ensembles and the frequency axis. And it will make the plots
% as described below. (Several are commented out.)
%
% This m-file operates on a matrix of ensembles of data in columns. You must specify the sampling rate that the
% data were obtained at as the second input to the function. It makes several plots where the first plot is the
% first two ensembles v. time. Then it looks at the histograms, auto correlations and power spectra for the first
% two ensembles, the average of the first ten ensembles and the mean of all of the available ensembles. It then
% plots the variance of the original signal v. the integral of the spectrum in frequency space and the
% autocorrelation at zero lag against the variance.
%
% (on hold since I will need to specify the x and y labels and other nicities. I suggest do this as post processing
% before you print it. You can specify the title for the plots as the third bit of information sent to this function.)
%
% USES: verticalfill.m
%
% NOTES: 
%		1. For slower computers, choose the length of the data record to be a power of two so the FFT is more efficient.
%		2. I am going to detrend the data in the ensembles first to get rid of any linear trend (single straight line fit).
%

[len,nE] = size(x);						% Get length of ensemble and number of ensembles.
t = [1:len]/sr; t=t';					% time axis
if exist('ndt')
	if lower(ndt) ~= 'n'
		x = detrend(x);
		disp('Linear trend were removed from the data')
	end
else		% Automatically detrend if ndt is not specified.
	x = detrend(x);
	disp('Linear trend were removed from the data')
end

% Default text strings
if ~exist('txxu')
	txxu = 'Hz'
end
if ~exist('txv')
	txv = 'x(t)';
end
if ~exist('txu')
	txu = 'units of x';
end
if ~exist('txx')
	txx = 'sample data';
end
txx = [', ' txx '.'];

% Set up some defaults for all figures.
scrsz=get(0,'ScreenSize'); % scrsz = [left bottom width height] in pixels
% Figures occupy a certain % of the screen and the full paper when printing
up = 50; rght = 0; phor = 1.00; pvert = 0.80;
set(0,'DefaultFigurePosition',[scrsz(1)+rght scrsz(2)+up phor*scrsz(3) pvert*scrsz(4)],'DefaultFigurePaperPosition',[0.25,0.25,8,10.5])
clear scrsz up rght phor pvert


% % Look at the first 8 ensembles.
% figure
% subplot(411),plot(t,x(:,1),'b',t,x(:,2),'r'); axis([0 max(t) -inf inf]); grid on; ylabel([txv ' [' txu ']'])
% title(['First 8 ' num2str(len) ' pt. ensembles of the ' txv ' data' txx]);
% subplot(412),plot(t,x(:,3),'b',t,x(:,4),'r'); axis([0 max(t) -inf inf]); grid on; ylabel([txv ' [' txu ']'])
% subplot(413),plot(t,x(:,5),'b',t,x(:,6),'r'); axis([0 max(t) -inf inf]); grid on; ylabel([txv ' [' txu ']'])
% subplot(414),plot(t,x(:,7),'b',t,x(:,8),'r'); axis([0 max(t) -inf inf]); grid on; ylabel([txv ' [' txu ']'])
% xlabel('t (sec)');

% % Compute the histograms using 20 bins NOTE that hist works on a matrix, psd does not.
% % Note, I have to take the mean of the histograms not the histogram of the mean of the signal.
% [H,c] = hist(x,50);						
% h10 = mean(H(:,[1:10]),2);		% Calculate the average of the first 10 ensembles
% h100=mean(H,2);								% Calculate the average of all the ensembles

% % Plot the signal and histograms
% figure
% subplot(4,1,1),plot(t,x(:,1),'b',t,x(:,2),'r'); xlabel('t (sec)'); axis([0 max(t) -inf inf]); grid on;
% title(['First 2 ' num2str(len) ' pt. ensembles of the ' txv ' data' txx]); ylabel([txv ' [' txu ']'])
% % Plot the Histogram
% subplot(4,1,2),plot(c,H(:,1),'bx',c,H(:,2),'xr',c,H(:,1),'b',c,H(:,2),'r'); grid on;
% title(' Histograms of 1^s^t 2 ensembles');
% subplot(4,1,3),plot(c,h10,'gx',c,h10,'g'); grid on;
% title('Mean of the histograms of the 1^s^t 10 ensembles');
% subplot(4,1,4),plot(c,h100,'gx',c,h100,'g'); grid on;
% % % Add a Gaussian curve - use the stats m-file that fits a scaled gaussian curve...
% % % Assume that randn gives a noise with std=1, which it seems to do.
% % hold on
% % G=1/sqrt(2*pi)*exp(-c.^2/2);
% % G=G*sqrt(2*pi)*nE;	% make the curve of simialr height to the samples by multiplying by the number of ensembles.
% % plot(c,G); hold off
% title(['Mean of the histograms of all ' num2str(nE) ' ensembles']); clear G;


% % Rxx - Calculate the autocovariance of the ensembles. Note that all lag vectors are the same.
% [R, lag] = xcorr(x(:,1));
% for i = 2:nE;
% 	rxx = xcorr(x(:,i));
% 	R=[R rxx];
% 	clear rxx
% end
% clear i
% lag = lag/sr;										% Make the lag axis relate to time by dividing by the sampling frequency.
% R = R/len;											% Normalizing the autocorrelation function by the number of samples, same as using 'biased'
% r10 = mean(R(:,[1:10]),2);			% Calculate the average of the first ten autocorrelations
% r100 = mean(R,2);								% Calculate the average of all the ensembles

% % Plot the Autocovariance/N (autocorrelation?)
% figure
% subplot(4,1,1),plot(t,x(:,1),'b',t,x(:,2),'r'); xlabel('t'); axis([0 max(t) -inf inf]); grid on;
% title(['First 2 ' num2str(len) ' pt. ensembles of the ' txv ' data' txx]); ylabel([txv ' [' txu ']'])
% subplot(4,1,2),plot(lag,R(:,1),'bx',lag,R(:,2),'r+'); axis([min(lag) max(lag) -inf inf]);
% title('Auto co-variance/N  of 1^s^t 2 ensembles'); grid on;
% subplot(4,1,3),plot(lag,r10,'g'); axis([min(lag) max(lag) -inf inf]);
% title('Mean of the Auto co-variance/N of the 1^s^t 10 ensembles'); grid on;
% ylabel(['Auto co-variance / N  [(' txu ')^2]   (Same label all 3 subplots)']);
% subplot(4,1,4),plot(lag,r100,'g'); axis([min(lag) max(lag) -inf inf]);
% title(['Mean of the Auto co-variance/N of all ' num2str(nE) ' ensembles']); grid on;
% xlabel('Lag (sec)')


% Sxx - Calculate the spectra for the segments: Note the frequency axes for all of the plots are the same.
% Since I am choosing weather or not to detrend above, do not detrend here again.
[S, F] = psd(x(:,1),len,sr,len,'none');
for i = 2:nE;
	s = psd(x(:,i),len,sr,len);
	S=[S s];
	clear s
end
% Nyquist frequency, in order to get the proper power spectra, the output from psd must be divided by the Nyquist
nyq=sr/2;
S=S/nyq;
clear i nyq
if nE >= 10
	s10 = mean(S(:,[1:10]),2);			% Calculate the average of the first ten spectra
end
s100 = mean(S,2);								% Calculate the average of all the ensembles

% ( From signal processing with Barry, SAR2.m and notes from his class.)
% Calculate the +/- 95% confidence bounds... create the curves in fig 3.10 of Jenkins and Watts
%		From the notes, the number of degrees of freedom (DOF) is the number of data segments * 2. (Zero overlap)
% The number of degrees of freedom is determined by the window (M = the width of the lag window) and the data length.
% here n is the length of the data segment we are interested in. Each piece of independent data has 2 DOF. If the
% data segments are overlapped, then there are formulae that can be found through Barry's notes.
% Old code
% Need to round it to the nearest integer to use chi2inv
%v=round(3*(len*nE)/len);
% New code
v=nE*2;	% DOF
% to be within a certain confidence bounds (1-a) mainly 0.95
a=0.05;
xer=[a/2 1-a/2];
% y(1) is upper bound on 95% confidence and y(2) is the lower one.
y=v./chi2inv(xer,v);
clear a xer
% Now add the 95% confidence interval to the mean line by multiplying the line I add by y(1) and plotting and by y(2)

% % Plot the spectra
% figure
% subplot(411),plot(t,x(:,1),'b',t,x(:,2),'r'); xlabel('t'); axis([0 max(t) -inf inf]); grid on;
% title(['First 2 ' num2str(len) ' pt. ensembles of the ' txv ' data' txx]); ylabel([txv ' [' txu ']'])
% subplot(412),semilogy(F,S(:,1),'b',F,S(:,2),'r'); axis([min(F) max(F) 0 inf]); grid on;
% title('PSD of 1^s^t 2 ensembles');
% subplot(413),semilogy(F,s10,'g'); axis([min(F) max(F) 0 inf]); grid on;
% title('Mean of the PSD''s of the 1^s^t 10 ensembles');
% ylabel(['Power Spectrum [(' txu ')^2 / Hz]   (same label for all 3 subplots)']);
% subplot(414),semilogy(F,s100,'g'); axis([min(F) max(F) 0 inf]); grid on;
% % The semilog plot said that there was negative Sxx data... I think it is in the axis specification... yup (fixed)
% % subplot(414),plot(F,s100,'g'); axis([min(F) max(F) -inf inf]); grid on; xlabel('Frequency [Hz]');
% xlabel('Frequency [Hz]');
% % try to estimate the noise level by taking the mean of the high frequency stuff and put it on the last plot.
% lin = mean(s100(F>=5)); ns = sqrt(lin*sr/len);
% hold on; plot( F(F>=5), lin*ones(size(F(F>=5))) ); hold off
% title(['Mean of the Power Spectra of all ' num2str(nE) ' ensembles. Noise est. = ' num2str(lin) ' or ' num2str(ns) ' [' txu ']']);
% clear lin


% Plot just the mean power spectrum with 95% confidence limits about it.
figure
if max(F) == 22050 % Option to make a special plot for spectra of the human audio range.
	% plot the data just for the audio range 20 - 20000 Hz
	ari = 1:length(F); ari = ari(F >= 20); ari = ari(F<=20000); %Audible range index
	% Omit 95% confidence limits for audio spectrum plots
	loglog(F(ari),s100(ari),'k');
	axis([min(F(ari)) max(F(ari)) 0 inf]);
else
	% Plot all frequencies with 95% confidence bounds
	semilogy(F,s100,'k', F,s100*y(1),'b.', F, s100*y(2),'b.');
%f	axis([min(F) max(F) 0 inf]);
 	axis([min(F) max(F) 0 inf]);
end
set(gca,'FontName','Times','FontSize',[16]); grid on;
ylabel(['Power Spectrum of ' txv ' [(' txu ')^2 / ' txxu ']']); xlabel(['Frequency [' txxu ']']);
title(['Mean of the Power Spectra of all ' num2str(nE) ' ensembles' txx]);
orient landscape


% % Summary plot for the average of all ensembles
% figure
% subplot(411),plot(t,x(:,1),'b',t,x(:,2),'r'); xlabel('t'); axis([0 max(t) -inf inf]); grid on;
% title(['First 2 ' num2str(len) ' pt. ensembles of the ' txv ' data' txx]); ylabel([txv ' [' txu ']'])
% subplot(412),plot(c,h100,'gx',c,h100,'g'); grid on;
% title(['Histogram of the mean of all ' num2str(nE) ' ensembles']);
% subplot(413),plot(lag,r100,'g'); axis([min(lag) max(lag) -inf inf]); xlabel('Lag (sec)'); grid on;
% title(['Mean of the Auto co-variance/N of all ' num2str(nE) ' ensembles']); ylabel(['Auto co-variance / N  [(' txu ')^2]']);
% subplot(414),semilogy(F,s100,'g'); axis([min(F) max(F) 0 inf]); grid on; xlabel('Frequency [Hz]');
% ylabel(['Power Spectrum [(' txu ')^2 / Hz]']);
% % Make a line on the plot with the mean of the spectrum from fc to the nyquist as an estimate of the noise level
% fc = 5;	% lower cuttoff frequency for the average
% lin = mean(s100(F>=fc)); ns = sqrt(lin*sr/len); % the average value and expressed as units of the signal
% hold on;	% 2 point line for the average.
% plot([min(F(F>=fc)) F(end)],[lin lin],'b',[min(F(F>=fc)) F(end)],[lin lin]*y(1),'b--',[min(F(F>=fc)) F(end)],[lin lin]*y(2),'b--');
% title(['Mean of the Power Spectra of all ' num2str(nE) ' ensembles. Mean for f > ' num2str(fc) ' [Hz] = ' num2str(lin) ' or ' num2str(ns) ' [' txu ']']);
% % Do a linear fit to the data above fc to see how far away it is from the mean.
% bb = regress(s100(F>=fc),[ones(size(F(F>=fc))) F(F>=fc)]);
% plot(F(F>=fc), bb(1)+bb(2)*F(F>=fc),'r');
% %plot(F(F>=fc),(bb(1)+bb(2)*F(F>=fc))*y(1),'r--', F(F>=fc),(bb(1)+bb(2)*F(F>=fc))*y(2),'r--'); % Add confidence limits
% % Maybe add a legend to reduce the data in the title? (not debugged)
% %s1 = ['S_x_x [' txu ']']; s2 = ['<S_x_x> > ' num2str(fc) ' [Hz]']; s3 = 'linear fit';
% %legend(s1,s2,0)
% hold off;
% clear lin ns


% % Plot the variance v. the integral of the spectrum
% figure
% V = var(x);
% iS = trapz(F,S);
% lin=[min(V) max(iS)];
% subplot(212), plot(V,iS,'.',lin,lin); grid on
% xlabel('variance of x (units of x)^2');
% ylabel('Integral of the spectrum of x');
% clear lin
% % Add the plot of the variance at zero lag
% Ro = max(R);
% lin=[min(V) max(Ro)];
% subplot(211), plot(V,Ro,'.',lin,lin); grid on;
% xlabel('variance of x (units of x)^2');
% ylabel('(Auto co-variance of x) / N at zero lag');
% title('Variance comparisons, Rxx(0) and the integral of the spectrum')
% clear lin
