Mercurial > hg > ltpda
diff testing/utp_1.1/utps/ao/utp_ao_lcpsd.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_lcpsd.m Tue Dec 06 18:42:11 2011 +0100 @@ -0,0 +1,842 @@ +% UTP_AO_LCPSD a set of UTPs for the ao/lcpsd method +% +% M Hewitson 06-08-08 +% +% $Id: utp_ao_lcpsd.m,v 1.29 2011/07/22 11:51:46 mauro Exp $ +% + +% <MethodDescription> +% +% The lcpsd method of the ao class computes the cross-spectral density between two +% time-series AOs on a log-frequency axis. +% +% </MethodDescription> + +function results = utp_ao_lcpsd(varargin) + + % Check the inputs + if nargin == 0 + + % Some keywords + class = 'ao'; + mthd = 'lcpsd'; + + 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_11(mthd, [at1 at1], ple1)]; % Test plotinfo doesn't disappear + results = [results utp_12]; % Test data lengths + results = [results utp_17]; % Test units handling: LCPSD + results = [results utp_18]; % Comparison with PSD + + 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 ~= 9, atest = false; end + % Check key + if ~io(3).plists.isparam('kdes'), atest = false; end + if ~io(3).plists.isparam('jdes'), atest = false; end + if ~io(3).plists.isparam('lmin'), 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('psll'), atest = false; end + if ~io(3).plists.isparam('times'), atest = false; end + if ~io(3).plists.isparam('split'), atest = false; end + % Check default value + if ~isequal(io(3).plists.find('kdes'), 100), atest = false; end + if ~isequal(io(3).plists.find('jdes'), 1000), atest = false; end + if ~isequal(io(3).plists.find('lmin'), 0), 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('psll'), 200), atest = false; end + if ~isEmptyDouble(io(3).plists.find('times')), atest = false; end + if ~isEmptyDouble(io(3).plists.find('split')), atest = false; end + % Check options + if ~isequal(io(3).plists.getOptionsForParam('kdes'), {100}), atest = false; end + if ~isequal(io(3).plists.getOptionsForParam('jdes'), {1000}), atest = false; end + if ~isequal(io(3).plists.getOptionsForParam('lmin'), {0}), atest = false; end + if ~isequal(io(3).plists.getOptionsForParam('win'), specwin.getTypes), atest = false; end + if ~isequal(io(3).plists.getOptionsForParam('olap'), {-1}), atest = false; end + if ~isequal(io(3).plists.getOptionsForParam('order'), {-1 0 1 2 3 4 5 6 7 8 9}), atest = false; end + if ~isequal(io(3).plists.getOptionsForParam('psll'), {200}), atest = false; end + if ~isequal(io(3).plists.getOptionsForParam('times'), {[]}), atest = false; end + if ~isequal(io(3).plists.getOptionsForParam('split'), {[]}), atest = false; 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 lcpsd 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 lcpsd method works for a vector of AOs as input. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + avec = [at1 at5]; + out = lcpsd(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 lcpsd method doesn't work with a matrix of AOs as input. + % + % </TestDescription> + function result = utp_03 + + % <SyntaxDescription> + % + % Test that the lcpsd method doesn't work for a matrix of AOs as input. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + amat = [at1 at5; at5 at6]; + out = lcpsd(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 lcpsd method works with a list of AOs as input. + % + % </TestDescription> + function result = utp_04 + + % <SyntaxDescription> + % + % Test that the lcpsd method works for a list of AOs as input. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + out = lcpsd(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 lcpsd method doesn't work with a mix of different shaped + % AOs as input. + % + % </TestDescription> + function result = utp_05 + + % <SyntaxDescription> + % + % Test that the lcpsd method doesn't work with an input of matrices and + % vectors and single AOs. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + out = lcpsd(at1,[at5 at6],at5,[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 lcpsd method properly applies history. + % + % </TestDescription> + function result = utp_06 + + % <SyntaxDescription> + % + % Test that the result of applying the lcpsd method can be processed back + % to an m-file. + % + % </SyntaxDescription> + + try + % <SyntaxCode> + out = lcpsd(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, 'lcpsd'), 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 lcpsd method can not modify the input AO. + % + % </TestDescription> + function result = utp_07 + + % <SyntaxDescription> + % + % Test that the lcpsd 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.lcpsd(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 lcpsd 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 = lcpsd(at5, at6); + out2 = lcpsd(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 lcpsd 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 lcpsd 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_12 + + % <TestDescription> + % + % Tests that differently sized data sets are treated properly + % + % </TestDescription> + function result = utp_12 + + % <SyntaxDescription> + % + % Test that applying lcpsd 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 lcpsd + win = 'Hanning'; + pl = plist('Win', win, 'order', -1); + out = lcpsd(a2,a1f,pl); + + stest = true; + % </SyntaxCode> + catch err + disp(err.message) + stest = false; + end + + % <AlgoDescription> + % + % 1) Check that lcpsd 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(1) ~= 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_12 + + %% 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) LCPSD 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) LCPSD 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 LCPSD 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) LCPSD of the white noise series + % + % Comparison with LPSD: + % 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) LCPSD of the white noise + % 6) LPSD 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 LCPSD(x,y) equals conj(LCPSD(y,x)) + % 2) Check that LCPSD(y,y) equals LPSD(y) + % 3) Check that LCPSD(x,x) equals LPSD(x) + + % + % </AlgoDescription> + + % <AlgoCode> + atest = true; + TOL = 1e-12; + + 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) + if any(abs(spec.Cxx.y - spec.S_1.y)./abs(spec.Cxx.y) > TOL) + atest = false; + end + end + if ne(spec.Cyy.data, spec.S_2.data) + if any(abs(spec.Cyy.y - spec.S_2.y)./abs(spec.Cyy.y) > TOL) + atest = false; + end + end + else + atest = false; + end + % </AlgoCode> + + % Return a result structure + result = utp_prepare_result(atest, stest, dbstack, mfilename); + end % END UTP_18 + + %% 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 lcpsd 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); + + otherwise + pl_spec = plist('Win', win_type, 'olap', olap, 'order', detrend_order); + + end + + if find(pli, 'win_obj') + % Calls the lcpsd applying the detrend and window internally + % (passig window object) + spec2.pl = pl_spec; + spec2.Cxy = lcpsd(a_1, a_2, spec2.pl); + spec2.Cyx = lcpsd(a_2, a_1, spec2.pl); + spec2.Cxx = lcpsd(a_1, a_1, spec2.pl); + spec2.Cyy = lcpsd(a_2, a_2, spec2.pl); + spec2.S_1 = simplifyYunits(lpsd(a_1, spec2.pl), ... + plist('prefixes', false, 'exceptions','Hz')); + spec2.S_2 = simplifyYunits(lpsd(a_2, spec2.pl), ... + plist('prefixes', false, 'exceptions','Hz')); + else + spec2 = struct; + end + % Calls the lcpsd applying the detrend and window internally + % (passig window name) + spec1.pl = pl_spec.pset('Win', win_type); + spec1.Cxy = lcpsd(a_1, a_2, spec1.pl); + spec1.Cyx = lcpsd(a_2, a_1, spec1.pl); + spec1.Cxx = lcpsd(a_1, a_1, spec1.pl); + spec1.Cyy = lcpsd(a_2, a_2, spec1.pl); + spec1.S_1 = simplifyYunits(lpsd(a_1, spec1.pl), ... + plist('prefixes', false, 'exceptions','Hz')); + spec1.S_2 = simplifyYunits(lpsd(a_2, spec1.pl), ... + plist('prefixes', false, 'exceptions','Hz')); + + end + +end