Mercurial > hg > ltpda
diff testing/utp_1.1/utps/ao/utp_ao_cpsd.m @ 44:409a22968d5e default
Add unit tests
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Tue, 06 Dec 2011 18:42:11 +0100 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/testing/utp_1.1/utps/ao/utp_ao_cpsd.m Tue Dec 06 18:42:11 2011 +0100 @@ -0,0 +1,1291 @@ +% UTP_AO_CPSD a set of UTPs for the ao/cpsd method +% +% M Hewitson 06-08-08 +% +% $Id: utp_ao_cpsd.m,v 1.43 2011/07/22 11:51:46 mauro Exp $ +% + +% <MethodDescription> +% +% The cpsd method of the ao class computes the cross-spectral density between two +% time-series AOs. +% +% </MethodDescription> + +function results = utp_ao_cpsd(varargin) + + % Check the inputs + if nargin == 0 + + % Some keywords + class = 'ao'; + mthd = 'cpsd'; + + results = []; + disp('******************************************************'); + disp(['**** Running UTPs for ' class '/' mthd]); + disp('******************************************************'); + + % Test AOs + [at1,at2,at3,at4,at5,at6] = eval(['get_test_objects_' class]); + + % Exception list for the UTPs: + [ple1,ple2,ple3,ple4,ple5,ple6] = get_test_ples(); + + % Get default window from the preferences + prefs = getappdata(0, 'LTPDApreferences'); + defaultWinType = char(prefs.getMiscPrefs.getDefaultWindow); + + % Run the tests + results = [results utp_01]; % getInfo call + results = [results utp_02]; % Vector input (only with two objects) + results = [results utp_03]; % Matrix input (not possible) + results = [results utp_04]; % List input (only with two objects) + results = [results utp_05]; % Test with mixed input (not possible) + results = [results utp_06]; % Test history is working + results = [results utp_07]; % Test the modify call works + results = [results utp_08]; % Test input data shape == output data shape + results = [results utp_09]; % Test output of the data + results = [results utp_10]; % Test against MATLAB's cpsd() + + results = [results utp_11(mthd, [at1 at1], ple1)]; % Test plotinfo doesn't disappear + + results = [results utp_17]; % Test units handling: CPSD + results = [results utp_18]; % Comparison with PSD + results = [results utp_24]; % Test data lengths + results = [results utp_25]; % Test Kaiser win and olap: CPSD + results = [results utp_51]; % Test number of averages: requested/obtained + results = [results utp_52]; % Test number of averages: correct number + results = [results utp_53]; % Test number of averages: syntax + + disp('Done.'); + disp('******************************************************'); + + elseif nargin == 1 % Check for UTP functions + if strcmp(varargin{1}, 'isutp') + results = 1; + else + results = 0; + end + else + error('### Incorrect inputs') + end + + %% UTP_01 + + % <TestDescription> + % + % Tests that the getInfo call works for this method. + % + % </TestDescription> + function result = utp_01 + + + % <SyntaxDescription> + % + % Test that the getInfo call works for no sets, all sets, and each set + % individually. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + % Call for no sets + io(1) = eval([class '.getInfo(''' mthd ''', ''None'')']); + % Call for all sets + io(2) = eval([class '.getInfo(''' mthd ''')']); + % Call for each set + for kk=1:numel(io(2).sets) + io(kk+2) = eval([class '.getInfo(''' mthd ''', ''' io(2).sets{kk} ''')']); + end + % </SyntaxCode> + stest = true; + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Check that getInfo call returned an minfo object in all cases. + % 2) Check that all plists have the correct parameters. + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % check we have minfo objects + if isa(io, 'minfo') + + % SET 'None' + if ~isempty(io(1).sets), atest = false; end + if ~isempty(io(1).plists), atest = false; end + % Check all Sets + if ~any(strcmpi(io(2).sets, 'Default')), atest = false; end + if numel(io(2).plists) ~= numel(io(2).sets), atest = false; end + % SET 'Default' + if io(3).plists.nparams ~= 8, atest = false; end + % Check key + if ~io(3).plists.isparam('nfft'), atest = false; end + if ~io(3).plists.isparam('win'), atest = false; end + if ~io(3).plists.isparam('olap'), atest = false; end + if ~io(3).plists.isparam('order'), atest = false; end + if ~io(3).plists.isparam('navs'), atest = false; end + if ~io(3).plists.isparam('times'), atest = false; end + if ~io(3).plists.isparam('split'), atest = false; end + if ~io(3).plists.isparam('psll'), atest = false; end + % Check default value + if ~isequal(io(3).plists.find('nfft'), -1), atest = false; end + if ~strcmpi(io(3).plists.find('win'), defaultWinType), atest = false; end + if ~isequal(io(3).plists.find('olap'), -1), atest = false; end + if ~isequal(io(3).plists.find('order'), 0), atest = false; end + if ~isequal(io(3).plists.find('navs'), -1), atest = false; end + if ~isEmptyDouble(io(3).plists.find('times')), atest = false; end + if ~isEmptyDouble(io(3).plists.find('split')), atest = false; end + if ~isequal(io(3).plists.find('psll'), 200), atest = false; end + % Check options + if ~isequal(io(3).plists.getOptionsForParam('nfft'), {-1}), atest = false; disp('1'); end + if ~isequal(io(3).plists.getOptionsForParam('win'), specwin.getTypes), atest = false;disp('2'); end + if ~isequal(io(3).plists.getOptionsForParam('olap'), {-1}), atest = false; disp('3');end + if ~isequal(io(3).plists.getOptionsForParam('order'), {-1 0 1 2 3 4 5 6 7 8 9}), atest = false;disp('4'); end + if ~isequal(io(3).plists.getOptionsForParam('navs'), {-1}), atest = false;disp('5'); end + if ~isequal(io(3).plists.getOptionsForParam('times'), {[]}), atest = false;disp('6'); end + if ~isequal(io(3).plists.getOptionsForParam('split'), {[]}), atest = false;disp('6'); end + if ~isequal(io(3).plists.getOptionsForParam('psll'), {200}), atest = false;disp('7'); end + end + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_01 + + %% UTP_02 + + % <TestDescription> + % + % Tests that the cpsd method works with a vector of AOs as input. (only + % with two objects in the vector) + % + % </TestDescription> + function result = utp_02 + + % <SyntaxDescription> + % + % Test that the cpsd method works for a vector of AOs as input. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + avec = [at1 at5]; + out = cpsd(avec); + % </SyntaxCode> + stest = true; + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Check that the number of elements in 'out' is equal to 1 + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % Check we have the correct number of outputs + if numel(out) ~= 1, atest = false; end + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_02 + + %% UTP_03 + + % <TestDescription> + % + % Tests that the cpsd method doesn't work with a matrix of AOs as input. + % + % </TestDescription> + function result = utp_03 + + % <SyntaxDescription> + % + % Test that the cpsd method doesn't work for a matrix of AOs as input. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + amat = [at1 at5 at6; at5 at6 at1]; + out = cpsd(amat); + % </SyntaxCode> + stest = false; + catch err + stest = true; + end + + % <AlgoDescription> + % + % 1) Nothing to check. + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_03 + + %% UTP_04 + + % <TestDescription> + % + % Tests that the cpsd method works with a list of AOs as input. + % + % </TestDescription> + function result = utp_04 + + % <SyntaxDescription> + % + % Test that the cpsd method works for a list of AOs as input. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + out = cpsd(at1,at5); + % </SyntaxCode> + stest = true; + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Check that the number of elements in 'out' is equal to 1 + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % Check we have the correct number of outputs + if numel(out) ~= 1, atest = false; end + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_04 + + %% UTP_05 + + % <TestDescription> + % + % Tests that the cpsd method doesn't work with a mix of different shaped + % AOs as input. + % + % </TestDescription> + function result = utp_05 + + % <SyntaxDescription> + % + % Test that the cpsd method doesn't work with an input of matrices and + % vectors and single AOs. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + out = cpsd([at5 at6],[at5 at1; at6 at1],at6); + stest = false; + % </SyntaxCode> + catch err + stest = true; + end + + % <AlgoDescription> + % + % 1) Nothing to check + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_05 + + %% UTP_06 + + % <TestDescription> + % + % Tests that the cpsd method properly applies history. + % + % </TestDescription> + function result = utp_06 + + % <SyntaxDescription> + % + % Test that the result of applying the cpsd method can be processed back + % to an m-file. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + out = cpsd(at5,at6); + mout = rebuild(out); + % </SyntaxCode> + stest = true; + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Check that the last entry in the history of 'out' corresponds to + % 'cpsd'. + % 2) Check that the re-built object is the same as 'out'. + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % Check the last step in the history of 'out' + if ~strcmp(out.hist.methodInfo.mname, 'cpsd'), atest = false; end + % Check the re-built object + if ~eq(mout, out, ple2), atest = false; end + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_06 + + %% UTP_07 + + % <TestDescription> + % + % Tests that the cpsd method can not modify the input AO. + % + % </TestDescription> + function result = utp_07 + + % <SyntaxDescription> + % + % Test that the cpsd method can not modify the input AO. + % The method must throw an error for the modifier call. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + % copy at1 to work with + ain = ao(at1); + % modify ain + ain.cpsd(at5); + % </SyntaxCode> + stest = false; + catch err + stest = true; + end + + % <AlgoDescription> + % + % 1) Nothing to check. + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_07 + + %% UTP_08 + + % <TestDescription> + % + % Test the shape of the output. + % + % </TestDescription> + function result = utp_08 + + % <SyntaxDescription> + % + % Test that the cpsd method keeps the data shape of the input object. The + % input AO must be an AO with row data and an AO with column data. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + out1 = cpsd(at5, at6); + out2 = cpsd(at6, at5); + % </SyntaxCode> + stest = true; + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Check that the shape of the output data doesn't change. + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % Check the shape of the output data + if size(out1.data.y, 2) ~= 1, atest = false; end + if size(out2.data.y, 1) ~= 1, atest = false; end + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_08 + + %% UTP_09 + + % <TestDescription> + % + % Check that the cpsd method pass back the output objects to a list of + % output variables or to a single variable. + % + % </TestDescription> + function result = utp_09 + + % <SyntaxDescription> + % + % This test is not longer necessary because the cpsd method pass back + % always only one object. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + % </SyntaxCode> + stest = true; + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Nothing to check. + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_09 + + %% UTP_10 + + % <TestDescription> + % + % Tests that the cpsd method agrees with MATLAB's cpsd when + % configured to use the same parameters. + % + % </TestDescription> + function result = utp_10 + + % <SyntaxDescription> + % + % Test that applying cpsd works on two AOs. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + % Construct two test AOs + nsecs = 10; + fs = 1000; + pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))'); + a1 = ao(pl); a2 = ao(pl); + % Filter one time-series + f2 = miir(plist('type', 'bandpass', 'fs', fs, 'order', 3, 'fc', [50 250])); + a1f = filter(a1, plist('filter', f2)); + % make some cross-power + a4 = a1f+a2; a4.setName; + % Compute cpsd + Nfft = 2*fs; + win = specwin('Hanning', Nfft); + pl = plist('Nfft', Nfft, 'Win', win.type, 'order', -1); + out = cpsd(a4,a1,pl); + % </SyntaxCode> + stest = true; + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Check that output agrees with the output of MATLAB's cpsd. + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % Compute cpsd using MATLAB's cpsd + [cxy, f] = cpsd(a4.data.y, a1.data.y, win.win, Nfft/2, Nfft, a1.data.fs); + if ~utils.math.isequal(cxy(:), out.data.y(:)) || ~utils.math.isequal(f, out.data.getX), atest = false; end + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_10 + + + %% UTP_17 + + % <TestDescription> + % + % Tests handling of units: + % 1) white noise produced from normal pdf, with a given mean value and + % sigma (distribution's 1st and 2nd orders) + % 2) white noise produced from normal pdf, with a given mean value and + % sigma (distribution's 1st and 2nd orders) + % 3) CPSD of the white noise series + % 4) compares the units of the input and output + % + + % </TestDescription> + function result = utp_17 + + % <SyntaxDescription> + % + % 1) Prepare the test tsdata: + % white noise from normal distribution + offset + % 2) Assign a random unit + % 3) Prepare the test tsdata: + % white noise from normal distribution + offset + % 4) Assign a random unit + % 5) CPSD of the white noise + % + % </SyntaxDescription> + + % <SyntaxCode> + try + + noise_type = 'Normal'; + win_type = 'BH92'; + + [a_1, a_2, spec, spec1] = prepare_analyze_noise(win_type, noise_type, plist); + + stest = true; + + catch err + disp(err.message) + stest = false; + end + % </SyntaxCode> + + % <AlgoDescription> + % + % 1) Check that (calculated CPSD yunits) equals + % input_1 units*input_2 units/Hz + % + % </AlgoDescription> + + % <AlgoCode> + atest = true; + u = simplifyYunits(a_1.* a_2, plist('prefixes', false, 'exceptions', 'Hz')); + if stest + if ne(spec.Cxy.yunits, u.yunits * unit('Hz^-1')) || ne(spec.Cxy.xunits, unit('Hz')) + atest = false; + end + if ne(spec.Cyx.yunits, u.yunits * unit('Hz^-1')) || ne(spec.Cyx.xunits, unit('Hz')) + atest = false; + end + else + atest = false; + end + % </AlgoCode> + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_17 + + %% UTP_18 + + % <TestDescription> + % + % Tests handling of units: + % 1) white noise produced from normal pdf, with a given mean value and + % sigma (distribution's 1st and 2nd orders) + % 2) white noise produced from normal pdf, with a given mean value and + % sigma (distribution's 1st and 2nd orders) + % 3) CPSD of the white noise series + % + % Comparison with PSD: + % 4) compares the off-diagonal terms to check they are complex-conjugated + % 5) compares the diagonal terms with PSD of the individual noise + % + + % </TestDescription> + function result = utp_18 + + % <SyntaxDescription> + % + % 1) Prepare the test tsdata: + % white noise from normal distribution + offset + % 2) Assign a random unit + % 3) Prepare the test tsdata: + % white noise from normal distribution + offset + % 4) Assign a random unit + % 5) CPSD of the white noise + % 6) PSD of the white noise + % + % </SyntaxDescription> + + % <SyntaxCode> + try + + noise_type = 'Uniform'; + win_type = 'BH92'; + + [a_1, a_2, spec, spec2] = prepare_analyze_noise(win_type, noise_type, plist); + + stest = true; + + catch err + disp(err.message) + stest = false; + end + % </SyntaxCode> + + % <AlgoDescription> + % + % 1) Check that CPSD(x,y) equals conj(CPSD(y,x)) + % 2) Check that CPSD(x,x) equals PSD(x) + % 3) Check that CPSD(y,y) equals PSD(y) + % + % </AlgoDescription> + + % <AlgoCode> + atest = true; + + if stest + if ne(spec.Cxy.y, conj(spec.Cyx.y)), atest = false; end + if ne(spec.Cxy.x, spec.Cyx.x), atest = false; end + if ne(spec.Cxx.data, spec.S_1.data), atest = false; end + if ne(spec.Cyy.data, spec.S_2.data), atest = false; end + else + atest = false; + end + % </AlgoCode> + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_18 + + + %% UTP_24 + + % <TestDescription> + % + % Tests that differently sized data sets are treated properly + % + % </TestDescription> + function result = utp_24 + + % <SyntaxDescription> + % + % Test that applying cpsd works on two AOs. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + % Construct two test AOs + nsecs = [10000:1:20000]; + fs = 1; + pl = plist('fs', fs, 'tsfcn', 'randn(size(t))'); + a1 = ao(pl.pset('nsecs', utils.math.randelement(nsecs, 1))); + a2 = ao(pl.pset('nsecs', utils.math.randelement(nsecs, 1))); + len_1 = a1.len; + len_2 = a2.len; + % Filter one time-series + f2 = miir(plist('type', 'bandpass', 'fs', fs, 'order', 3, 'fc', [.050 .25])); + a1f = filter(a1, plist('filter', f2)); + % Compute cpsd + Nfft = -1; + win = 'Hanning'; + pl = plist('Nfft', Nfft, 'Win', win, 'order', -1); + out = cpsd(a2,a1f,pl); + % </SyntaxCode> + stest = true; + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Check that cpsd used the length of the shortest ao. + % + % </AlgoDescription> + + atest = true; + if stest + % <AlgoCode> + % Compare the nfft with the length of the input data + + if out.x(2) ~= 1/min(len_1,len_2) + atest = false; + end + % </AlgoCode> + else + atest = false; + end + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_24 + + %% UTP_25 + + % <TestDescription> + % + % Tests handling of units: + % 1) white noise produced from normal pdf, with a given mean value and + % sigma (distribution's 1st and 2nd orders) + % 2) white noise produced from normal pdf, with a given mean value and + % sigma (distribution's 1st and 2nd orders) + % 3) CPSD of the white noise series + % 4) compares the units of the input and output + % + + % </TestDescription> + function result = utp_25 + + % <SyntaxDescription> + % + % 1) Prepare the test tsdata: + % white noise from normal distribution + offset + % 2) Assign a random unit + % 3) Prepare the test tsdata: + % white noise from normal distribution + offset + % 4) Assign a random unit + % 5) CPSD of the white noise + % + % </SyntaxDescription> + + % <SyntaxCode> + try + + % Build time-series test data + fs = 1; + nsecs = 86400; + sigma_distr_1 = 4.69e-12; + mu_distr_1 = -5.11e-14; + sigma_distr_2 = 6.04e-9; + mu_distr_2 = 1.5e-10; + + % White noise + type = 'Normal'; + + a_n = ao(plist('waveform', 'noise', ... + 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1)); + a_const = ao(mu_distr_1); + a_1 = a_n + a_const; + + a_n = ao(plist('waveform', 'noise', ... + 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_2)); + a_const = ao(mu_distr_2); + a_2 = a_n + a_const; + + % Set units and prefix from those supported + unit_list = unit.supportedUnits; + % remove the first empty unit '' from the list, because then is it + % possible that we add a prefix to an empty unit + unit_list = unit_list(2:end); + prefix_list = unit.supportedPrefixes; + a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))])); + a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))])); + + % Evaluate the cpsd of the time-series data, using Kaiser window + % Psll and olap are not set + win = ('Kaiser'); + detrend = 0; + n_pts = nsecs*fs/10; + + C = cpsd(a_1, a_2, plist('Win', win, 'Nfft', n_pts, 'order', detrend)); + + stest = true; + + catch err + disp(err.message) + stest = false; + end + % </SyntaxCode> + + % <AlgoDescription> + % + % 1) Check that (calculated CPSD yunits) equals + %input_1 units*input_2 units/Hz + + % </AlgoDescription> + + % <AlgoCode> + atest = true; + u = simplifyYunits(a_1.* a_2, plist('prefixes', false, 'exceptions', 'Hz')); + if stest + if ne(C.yunits, u.yunits * unit('Hz^-1')) || ne(C.xunits, unit('Hz')) + atest = false; + end + else + atest = false; + end + % </AlgoCode> + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_25 + + %% UTP_51 + + % <TestDescription> + % + % Tests the possibility to set the number of averages rather than setting the Nfft: + % 1) white noise produced from normal pdf, with: + % a given mean value and sigma (distribution's 1st and 2nd order) + % 2) cpsd of the noise, without detrending, random window, set number of + % averages + % 3) check the effective number of averages + % + + % </TestDescription> + function result = utp_51 + + % <SyntaxDescription> + % + % 1) Prepare the test tsdata: + % white noise from normal distribution + offset + % 2) cpsd of the noise, without detrending, random window, set number of + % averages + % + % </SyntaxDescription> + + % <SyntaxCode> + try + + noise_type = 'Normal'; + + % Evaluate the cpsd of the white noise time-series data + [a_1, a_2, C1, C2, navs] = prepare_analyze_noise_navs(noise_type, plist); + + stest = true; + + catch err + disp(err.message) + stest = false; + end + % </SyntaxCode> + + % <AlgoDescription> + % + % 1) Check that calculated navs are identical to those requested + % + % </AlgoDescription> + + % <AlgoCode> + atest = true; + + if stest + if ne(navs, C1.data.navs) + if ne(find(C1.hist.plistUsed, 'navs'), C1.data.navs) + atest = false; + end + end + else + atest = false; + end + % </AlgoCode> + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_51 + + %% UTP_52 + + % <TestDescription> + % + % Tests the possibility to set the number of averages rather than setting the Nfft: + % 1) white noise produced from uniform pdf, with: + % a given mean value and sigma (distribution's 1st and 2nd order) + % 2) cpsd of the noise, without detrending, random window, random navs + % 3) get the number of averages + % 4) get the nfft used + % 5) run cpsd again, with the nfft used + % 6) compare the calculated objects + % + + % </TestDescription> + function result = utp_52 + + % <SyntaxDescription> + % + % 1) white noise produced from uniform pdf, with: + % a given mean value and sigma (distribution's 1st and 2nd order) + % 2) cpsd of the noise, without detrending, random window, random navs + % 3) get the number of averages + % 4) get the nfft used + % 5) run cpsd again, with the nfft used + % + % </SyntaxDescription> + + % <SyntaxCode> + try + + noise_type = 'Uniform'; + + % Evaluate the cpsd of the white noise time-series data + [a_1, a_2, C1, C2, navs] = prepare_analyze_noise_navs(noise_type, plist); + + stest = true; + + catch err + disp(err.message) + stest = false; + end + % </SyntaxCode> + + % <AlgoDescription> + % + % 1) Check that calculated objects C1 and C2 are identical + % + % </AlgoDescription> + + % <AlgoCode> + atest = true; + + if stest + % Compare the output objects + if ne(C1, C2, ple3) + atest = false; + end + else + atest = false; + end + % </AlgoCode> + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_52 + + %% UTP_53 + + % <TestDescription> + % + % Tests the possibility to set the number of averages rather than setting the Nfft: + % 1) white noise produced from normal pdf, with: + % a given mean value and sigma (distribution's 1st and 2nd order) + % 2) cpsd of the noise, without detrending, random window, random navs + % 3) get the number of averages + % 4) get the nfft used + % 5) run cpsd again, with the nfft used + % 6) compare navs, nfft, psds + % + + % </TestDescription> + function result = utp_53 + + % <SyntaxDescription> + % + % 1) white noise produced from normal pdf, with: + % a given mean value and sigma (distribution's 1st and 2nd order) + % 2) cpsd of the noise, without detrending, random window, random navs + % 3) get the number of averages + % 4) get the nfft used + % 5) run cpsd again, with the nfft used + % 6) run cpsd again, with conflicting parameters, and verify it uses + % nfft rather than navs + % + % </SyntaxDescription> + + % <SyntaxCode> + try + + noise_type = 'Uniform'; + + % Evaluate the cpsd of the white noise time-series data + [a_1, a_2, C1, C2, navs] = prepare_analyze_noise_navs(noise_type, plist); + + npts_3 = fix(find(C1.hist.plistUsed, 'Nfft')/2); + + % Calculates the cpsd asking for the number of points AND the window length + pl_spec = C1.hist.plistUsed; + pl_spec.pset('Nfft', npts_3, 'navs', navs); + C3 = cpsd(a_1, a_2, pl_spec); + + stest = true; + + catch err + disp(err.message) + stest = false; + end + % </SyntaxCode> + + % <AlgoDescription> + % + % 1) Check that calculated objects C1 and C2 are identical + % 2) Check that C3 used different values + % + % </AlgoDescription> + + % <AlgoCode> + atest = true; + + if stest + % Compare the navs written in the output object with the requested one + if ne(C1,C2,ple3) || ... + ne(find(C3.hist.plistUsed, 'Nfft'), npts_3) || eq(C3.data.navs, navs) + atest = false; + end + else + atest = false; + end + % </AlgoCode> + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_53 + + %% Helper function for window call construction + + function [a_1, a_2, spec1, spec2] = prepare_analyze_noise(win_type, noise_type, pli) + % Array of parameters to pick from + fs_list = [0.1;1;2;5;10]; + nsecs_list = [20 100 1000:1000:10000]'; + sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]'; + trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]'; + + + % Build time-series test data + + % Picks the values at random from the list + fs = utils.math.randelement(fs_list, 1); + nsecs = utils.math.randelement(nsecs_list, 1); + sigma_distr_1 = utils.math.randelement(sigma_distr_list, 1); + sigma_distr_2 = utils.math.randelement(sigma_distr_list, 1); + trend_0_1 = utils.math.randelement(trend_0_list, 1); + trend_0_2 = utils.math.randelement(trend_0_list, 1); + + % Pick units and prefix from those supported + unit_list = unit.supportedUnits; + % remove the first empty unit '' from the list, because then is it + % possible that we add a prefix to an empty unit + unit_list = unit_list(2:end); + prefix_list = unit.supportedPrefixes; + + % White noise + a_n = ao(plist('waveform', 'noise', ... + 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1)); + + % Constant signal + a_c = ao(trend_0_1); + + % Total signal + a_1 = a_n + a_c; + + % White noise + a_n = ao(plist('waveform', 'noise', ... + 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_2)); + % Constant signal + a_c = ao(trend_0_2); + + % Total signal + a_2 = a_n + a_c; + + % Set units + a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))])); + a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))])); + + % Evaluate the cpsd of the white noise time-series data + olap = 0; + detrend_order = 0; + + switch lower(win_type) + case 'kaiser' + psll = find(pli, 'psll'); + if isempty(psll) + psll = find(ao.getInfo('psd').plists, 'psll'); + end + pl_spec = plist('Win', win_type, 'psll', psll, 'olap', olap, 'order', detrend_order); + + case 'levelledhanning' + levelCoef = find(pli, 'levelCoef'); + if isempty(levelCoef) + levelCoef = 1; + end + pl_spec = plist('Win', win_type, 'levelCoef', levelCoef, 'olap', olap, 'order', detrend_order); + + otherwise + pl_spec = plist('Win', win_type, 'olap', olap, 'order', detrend_order); + + end + + if find(pli, 'win_obj') + % Calls the cpsd applying the detrend and window internally + % (passig window object) + spec2.pl = pl_spec; + spec2.Cxy = cpsd(a_1, a_2, spec2.pl); + spec2.Cyx = cpsd(a_2, a_1, spec2.pl); + spec2.Cxx = cpsd(a_1, a_1, spec2.pl); + spec2.Cyy = cpsd(a_2, a_2, spec2.pl); + spec2.S_1 = simplifyYunits(psd(a_1, spec2.pl), ... + plist('prefixes', false, 'exceptions','Hz')); + spec2.S_2 = simplifyYunits(psd(a_2, spec2.pl), ... + plist('prefixes', false, 'exceptions','Hz')); + else + spec2 = struct; + end + % Calls the cpsd applying the detrend and window internally + % (passig window name) + spec1.pl = pl_spec.pset('Win', win_type); + spec1.Cxy = cpsd(a_1, a_2, spec1.pl); + spec1.Cyx = cpsd(a_2, a_1, spec1.pl); + spec1.Cxx = cpsd(a_1, a_1, spec1.pl); + spec1.Cyy = cpsd(a_2, a_2, spec1.pl); + spec1.S_1 = simplifyYunits(psd(a_1, spec1.pl), ... + plist('prefixes', false, 'exceptions','Hz')); + spec1.S_2 = simplifyYunits(psd(a_2, spec1.pl), ... + plist('prefixes', false, 'exceptions','Hz')); + + end + + %% Helper function for window call construction, navs option + + function [a_1, a_2, C1, C2, navs] = prepare_analyze_noise_navs(noise_type, pli) + % Array of parameters to pick from + fs_list = [0.1;1;2;5;10]; + nsecs_list = [2000:1000:10000]'; + sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]'; + trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]'; + + % Build time-series test data + + % Picks the values at random from the list + fs = utils.math.randelement(fs_list, 1); + nsecs = utils.math.randelement(nsecs_list, 1); + sigma_distr = utils.math.randelement(sigma_distr_list, 1); + trend_0 = utils.math.randelement(trend_0_list, 1); + + % Pick units and prefix from those supported + unit_list = unit.supportedUnits; + % remove the first empty unit '' from the list, because then is it + % possible that we add a prefix to an empty unit + unit_list = unit_list(2:end); + prefix_list = unit.supportedPrefixes; + + % White noise + a_n1 = ao(plist('waveform', 'noise', ... + 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr)); + a_n2 = ao(plist('waveform', 'noise', ... + 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr)); + + % Constant signal + a_c = ao(trend_0); + + % Total signals + a_1 = a_n1 + a_c; + a_2 = a_n2 + a_c; + + % Set units + a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))])); + a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))])); + + % Evaluate the cpsd of the white noise time-series data + olap = 0; + detrend_order = 0; + n_pts = -1; + + navs = fix(utils.math.randelement(logspace(0,log10(max(0,a_1.len/10)),50),1)); + + % Evaluate the cpsd of the white noise time-series data + % Window + win_list = specwin.getTypes; + win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1); + win_type = win_type{1}; + + switch lower(win_type) + case 'kaiser' + psll = utils.math.randelement([0:10:200],1); + if psll == 0 + psll = find(ao.getInfo('psd').plists, 'psll'); + end + pl_spec = plist('Win', win_type, 'psll', psll, 'olap', olap, 'order', detrend_order); + otherwise + pl_spec = plist('Win', win_type, 'olap', olap, 'order', detrend_order); + end + + % Calls cpsd asking for the number of averages + pl_spec.pset('Nfft', n_pts, 'navs', navs); + C1 = cpsd(a_1, a_2, pl_spec); + + % Calls cpsd asking for the number of points just evaluated + pl_spec.pset('Nfft', find(C1.hist.plistUsed, 'Nfft')); + pl_spec.remove('navs'); + C2 = cpsd(a_1, a_2, pl_spec); + + end + +end