% multispectest.m
%
% 13 Jan 2004 (David S. Ciochetto)
% Test multispec.m
%
% Routine to generate a fake signal to study power spectra estiamtes in multispec.m
% Generates nE ensembles of the same deterministic signal with gaussian white noise.
% To specify a variance uncomment the randvec.m and xx parts.
%
% USES: randvec.m
% 
%clear all; close all;

%%% INPUT %%%
nE  = 10;						% Number of ensembles
len = 2^9;					% Length of the data record (samples)
sr  = 64;						% Sampling Rate (e.g. 64 Hz = 64 samples/s)
su  = 'cycles/day';					% Units on sampling rate
yu  = 'm';					% units of the fake signals
amp = 1;						% Main Amplitude (units of x)
np  = 1.0;					% Percentage of noise with respect to the signal amplitude as a decimal
detrnd = 'y';				% Have multispec remove a linear trend? ('y'/'n')

	% Deterministic signal (2 sinusoidal components only unless changed below.) No trend.
		determ = 'y';		% Add the deterministic part? ('y'/'n')
		% Sinusoid(s):
		sf1 = 0.47;	 		% Frequency (Hz)
		sf2 = 2.65;
		ap1 = 1*amp;		% Amplitude of signal 1
		ap2 = 0.5*amp;	% Amplitude of signal 2


%%%%%%%% no more user input below this %%%%%%%%%
x = randn(len,nE);						% generate nE random white noise time series of length len
t = [1:len]/sr; t=t';					% time axis
noise = np*amp*x;							% Adjust the magnitude of the noise here
x = noise;
sig1 = ap1*sin(sf1*2*pi*t);		% Generate the deterministic signal
sig2 = ap2*sin(sf2*2*pi*t);
if determ == 'y'							% Add the deterministic signal if desired.
	sig = sig1 + sig2;
else
	sig = zeros(len,1);
end
for i = 1:nE									% Construct the signal
	x(:,i) = sig + noise(:,i);
end

% Nice text descriptions of the fake signal
if determ == 'y'
	%ytxt = [num2str(ap1) 'Sin(' num2str(sf1*2*pi) 't) + ' num2str(ap2) 'Sin(' num2str(sf2*2*pi) 't) + noise'];
	ytxt = 'MLD';
else
	ytxt = 'white noise'; 
end
clear i np amp ap1 ap2 sf1 sf2 sig sig1 sig2 noise determ

vx = mean(var(x));
multispec(x,sr,detrnd,su,ytxt,yu,['Test of spectral estimator, var of input = ' num2str(vx)]);

% % another plot with noise generated with randvec to specify the variance.
% xx=randvec(nE,ones(len,1)*0.4,0.005)';
% vxx = mean(var(xx));
% multispec(xx,sr,'y','white noise','m/s',['Test of spectral estimator, var of input = ' num2str(vxx)]);
