function varargout = pwelch(sys,varargin)
% ds = pwelch(sys,...)
%
% Simulates spectrum of sys using pwelch. If called without output
% arguments the resulting spectra are ploted. The output ds is a dadac
% object containing the simulated (colored) spectra and the white
% noise input (measurement noise). If no parameter structure is given,
% model sample frequency must be specified via optional argument.
%
% Optional arguments (key/value pairs)
% - 'reference': show SN spectrum
% - 'plots': always show plots
% Note that optional arguments have lower precedence than parameter
% structure!

% TODO:
%  - change postprocessing!!
%  - should this automatically postprocess?

% NOTE:
%  - this uses dadac objects!

% PSD can only be calculated for noisy systems
if ~sys.isnoisy()
    error('system must be noisy')
end

%% instance dac object
dao = dac.instance();

%% parse input arguments
args = varargin;
nargs = length(args);

% defaults for optional arguments
PLOT_REFERENCE_SPECTRUM = false;
SHOW_PLOTS = false;

% parse rest of varargs
% arguments from pars structure take precedence over these
for ii = 1:2:nargs
    switch args{ii}
      case 'plots', SHOW_PLOTS = args{ii+1};
      case 'reference', PLOT_REFERENCE_SPECTRUM = args{ii+1};
      otherwise, warning(['unrecognized argument ' args{ii}])
    end
end

%% check input

% check if Ntsim is defined
if isempty(dao.Ntsim)
    error('Ntsim must be set')
end


%% local configuration
if ~dao.quiet ddisp = @disp; else ddisp = @() disp(''); end

% initialize PRNG
if dao.prngseed
    rand('state', dao.prngseed);
    randn('state', dao.prngseed);
end

%% main
if sys.Ts == 0
    dt = 1/dao.Fs;
else
    dt = sys.Ts;
end

% sample trajectory
xinit = zeros(size(sys.a,1),1);
tvec = (0:dt:(dao.Ntsim-1)*dt)';

% sample measurement + white noise vacuum spectrum (= obsnoise)
simdata = dadac(2,'type','ts');
[simdata(1).yData,~,~,~,simdata(2).yData] = initial(sys,xinit,tvec);
simdata(1).label = 'simulated measurement';
simdata(2).label = 'simulated measurement white noise';

% postprocess white noise if required
% TODO: CHANGE THIS
if  dao.postprocess && ~strcmp(dao.postp.method,'none')
    bpm = bandpass_model(dao.postp.fdesign,dao.postp.method);
    bpm = append(bpm,bpm);
    simdata(2).yData = lsim(ss(bpm),simdata(2).yData,tvec);
end

% calculate spectrum
frdata = pwelch(simdata);

%% asign output
% plot spectra
for ii = 1:numel(frdata)
    if nargout == 0 || SHOW_PLOTS
        pwelchplot(frdata(ii));
    else
        varargout{1} = frdata;
    end
end