Mercurial > hg > ltpda

comparison testing/utp_1.1/utps/ao/utp_ao_cohere.m @ 44:409a22968d5e default

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add unit tests
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Tue, 06 Dec 2011 18:42:11 +0100
parents
children
comparison
equal deleted inserted replaced
43:bc767aaa99a8 44:409a22968d5e
1 % UTP_AO_COHERE a set of UTPs for the ao/cohere method
2 %
3 % M Hewitson 06-08-08
4 %
5 % $Id: utp_ao_cohere.m,v 1.44 2011/07/22 12:29:58 mauro Exp $
6 %
7
8 % <MethodDescription>
9 %
10 % The cohere method of the ao class computes the coherence between two
11 % time-series AOs.
12 %
13 % </MethodDescription>
14
15 function results = utp_ao_cohere(varargin)
16
17 % Check the inputs
18 if nargin == 0
19
20 % Some keywords
21 class = 'ao';
22 mthd = 'cohere';
23
24 results = [];
25 disp('******************************************************');
26 disp(['**** Running UTPs for ' class '/' mthd]);
27 disp('******************************************************');
28
29 % Test AOs
30 [at1,at2,at3,at4,at5,at6] = eval(['get_test_objects_' class]);
31
32 % Exception list for the UTPs:
33 [ple1,ple2,ple3,ple4,ple5,ple6] = get_test_ples();
34
35 % Get default window from the preferences
36 prefs = getappdata(0, 'LTPDApreferences');
37 defaultWinType = char(prefs.getMiscPrefs.getDefaultWindow);
38
39 % Run the tests
40 results = [results utp_01]; % getInfo call
41 results = [results utp_02]; % Vector input (only with two objects)
42 results = [results utp_03]; % Matrix input (not possible)
43 results = [results utp_04]; % List input (only with two objects)
44 results = [results utp_05]; % Test with mixed input (not possible)
45 results = [results utp_06]; % Test history is working
46 results = [results utp_07]; % Test the modify call works
47 results = [results utp_08]; % Test input data shape == output data shape
48 results = [results utp_09]; % Test output of the data
49 results = [results utp_10]; % Test the basic usage against MATLAB mscohere
50
51 results = [results utp_11(mthd, [at1 at1], ple1)]; % Test plotinfo doesn't disappear
52
53 results = [results utp_12]; % Test basic symmetry properties of cohere (C)
54 results = [results utp_13]; % Test basic symmetry properties of cohere (MS)
55 results = [results utp_14]; % Test basic symmetry properties of cohere (C)
56 results = [results utp_15]; % Test basic symmetry properties of cohere (MS)
57 results = [results utp_16]; % Test basic relationship (MS) <-> (C)
58 results = [results utp_17]; % Test units handling: complex cohere
59 results = [results utp_18]; % Test units handling: magnitude-squared cohere
60 results = [results utp_19]; % Test data lengths
61 results = [results utp_20]; % Test with single window
62 results = [results utp_21]; % Test number of averages: requested/obtained
63 results = [results utp_22]; % Test number of averages: correct number
64 results = [results utp_23]; % Test number of averages: syntax
65 results = [results utp_24]; % Test the basic usage against MATLAB mscohere
66 results = [results utp_25]; % Test Kaiser win and olap: (C)
67 results = [results utp_26]; % Test Kaiser win and olap: (MS)
68 results = [results utp_30]; % Special cases: same input
69
70 disp('Done.');
71 disp('******************************************************');
72
73 elseif nargin == 1 % Check for UTP functions
74 if strcmp(varargin{1}, 'isutp')
75 results = 1;
76 else
77 results = 0;
78 end
79 else
80 error('### Incorrect inputs')
81 end
82
83 %% UTP_01
84
85 % <TestDescription>
86 %
87 % Tests that the getInfo call works for this method.
88 %
89 % </TestDescription>
90 function result = utp_01
91
92
93 % <SyntaxDescription>
94 %
95 % Test that the getInfo call works for no sets, all sets, and each set
96 % individually.
97 %
98 % </SyntaxDescription>
99
100 try
101 % <SyntaxCode>
102 % Call for no sets
103 io(1) = eval([class '.getInfo(''' mthd ''', ''None'')']);
104 % Call for all sets
105 io(2) = eval([class '.getInfo(''' mthd ''')']);
106 % Call for each set
107 for kk=1:numel(io(2).sets)
108 io(kk+2) = eval([class '.getInfo(''' mthd ''', ''' io(2).sets{kk} ''')']);
109 end
110 % </SyntaxCode>
111 stest = true;
112 catch err
113 disp(err.message)
114 stest = false;
115 end
116
117 % <AlgoDescription>
118 %
119 % 1) Check that getInfo call returned an minfo object in all cases.
120 % 2) Check that all plists have the correct parameters.
121 %
122 % </AlgoDescription>
123
124 atest = true;
125 if stest
126 % <AlgoCode>
127 % check we have minfo objects
128 if isa(io, 'minfo')
129
130 % SET 'None'
131 if ~isempty(io(1).sets), atest = false; end
132 if ~isempty(io(1).plists), atest = false; end
133 % Check all Sets
134 if ~any(strcmpi(io(2).sets, 'Default')), atest = false; end
135 if numel(io(2).plists) ~= numel(io(2).sets), atest = false; end
136 % SET 'Default'
137 if io(3).plists.nparams ~= 9, atest = false; end
138 % Check key
139 if ~io(3).plists.isparam('nfft'), atest = false; end
140 if ~io(3).plists.isparam('win'), atest = false; end
141 if ~io(3).plists.isparam('olap'), atest = false; end
142 if ~io(3).plists.isparam('type'), atest = false; end
143 if ~io(3).plists.isparam('order'), atest = false; end
144 if ~io(3).plists.isparam('navs'), atest = false; end
145 if ~io(3).plists.isparam('times'), atest = false; end
146 if ~io(3).plists.isparam('split'), atest = false; end
147 if ~io(3).plists.isparam('psll'), atest = false; end
148 % Check default value
149 if ~isequal(io(3).plists.find('nfft'), -1), atest = false; end
150 if ~strcmpi(io(3).plists.find('win'), defaultWinType), atest = false; end
151 if ~isequal(io(3).plists.find('olap'), -1), atest = false; end
152 if ~isequal(io(3).plists.find('type'), 'C'), atest = false; end
153 if ~isequal(io(3).plists.find('order'), 0), atest = false; end
154 if ~isequal(io(3).plists.find('navs'), -1), atest = false; end
155 if ~isEmptyDouble(io(3).plists.find('times')), atest = false; end
156 if ~isEmptyDouble(io(3).plists.find('split')), atest = false; end
157 if ~isequal(io(3).plists.find('psll'), 200), atest = false; end
158 % Check options
159 if ~isequal(io(3).plists.getOptionsForParam('nfft'), {-1}), atest = false; end
160 if ~isequal(io(3).plists.getOptionsForParam('win'), specwin.getTypes), atest = false; end
161 if ~isequal(io(3).plists.getOptionsForParam('olap'), {-1}), atest = false; end
162 if ~isequal(io(3).plists.getOptionsForParam('type'), {'C', 'MS'}), atest = false; end
163 if ~isequal(io(3).plists.getOptionsForParam('order'), {-1 0 1 2 3 4 5 6 7 8 9}), atest = false; end
164 if ~isequal(io(3).plists.getOptionsForParam('navs'), {-1}), atest = false; end
165 if ~isequal(io(3).plists.getOptionsForParam('times'), {[]}), atest = false; end
166 if ~isequal(io(3).plists.getOptionsForParam('split'), {[]}), atest = false; end
167 if ~isequal(io(3).plists.getOptionsForParam('psll'), {200}), atest = false; end
168 end
169 % </AlgoCode>
170 else
171 atest = false;
172 end
173
174 % Return a result structure
175 result = utp_prepare_result(atest, stest, dbstack, mfilename);
176 end % END UTP_01
177
178 %% UTP_02
179
180 % <TestDescription>
181 %
182 % Tests that the cohere method works with a vector of AOs as input. (only
183 % with two objects in the vector)
184 %
185 % </TestDescription>
186 function result = utp_02
187
188 % <SyntaxDescription>
189 %
190 % Test that the cohere method works for a vector of AOs as input.
191 %
192 % </SyntaxDescription>
193
194 try
195 % <SyntaxCode>
196 avec = [at1 at5];
197 out = cohere(avec);
198 % </SyntaxCode>
199 stest = true;
200 catch err
201 disp(err.message)
202 stest = false;
203 end
204
205 % <AlgoDescription>
206 %
207 % 1) Check that the number of elements in 'out' is equal to 1.
208 % 2) Check that each output AO contains the correct data.
209 %
210 % </AlgoDescription>
211
212 atest = true;
213 if stest
214 % <AlgoCode>
215 % Check we have the correct number of outputs
216 if numel(out) ~= 1, atest = false; end
217
218 TOL = 1e-13;
219
220 % Get shortest vector
221 lmin = min([length(at1.y), length(at5.y), length(at6.y)]);
222 % Set Nfft
223 Nfft = lmin;
224 % Get default window
225 if strcmpi(defaultWinType, 'kaiser')
226 win = specwin(defaultWinType, Nfft, find(ao.getInfo('cohere').plists, 'psll'));
227 else
228 win = specwin(defaultWinType, Nfft);
229 end
230 % Compute magnitude squared coherence estimate with MATLAB
231 % out: at1->at5
232 [cxy, f] = mscohere(at1.y(1:lmin), at5.y(1:lmin), win.win, Nfft/2, Nfft, at1.fs);
233 if any(abs(out.y-cxy > TOL)), atest = false; end
234 if any(abs(out.x-f > TOL)), atest = false; end
235 % </AlgoCode>
236 else
237 atest = false;
238 end
239
240 % Return a result structure
241 result = utp_prepare_result(atest, stest, dbstack, mfilename);
242 end % END UTP_02
243
244 %% UTP_03
245
246 % <TestDescription>
247 %
248 % Test that the cohere method doesn't work for a matrix of AOs as input.
249 %
250 % </TestDescription>
251 function result = utp_03
252
253 % <SyntaxDescription>
254 %
255 % Test that the cohere method doesn't work for a matrix of AOs as input.
256 %
257 % </SyntaxDescription>
258
259 try
260 % <SyntaxCode>
261 amat = [at1 at2;at5 at6];
262 out = cohere(amat);
263 % </SyntaxCode>
264 stest = false;
265 catch err
266 stest = true;
267 end
268
269 % <AlgoDescription>
270 %
271 % 1) Nothing to check.
272 %
273 % </AlgoDescription>
274
275 atest = true;
276 if stest
277 % <AlgoCode>
278 % </AlgoCode>
279 else
280 atest = false;
281 end
282
283 % Return a result structure
284 result = utp_prepare_result(atest, stest, dbstack, mfilename);
285 end % END UTP_03
286
287 %% UTP_04
288
289 % <TestDescription>
290 %
291 % Tests that the cohere method works with a list of AOs as input.
292 %
293 % </TestDescription>
294 function result = utp_04
295
296 % <SyntaxDescription>
297 %
298 % Test that the cohere method works for a list of AOs as input.
299 %
300 % </SyntaxDescription>
301
302 try
303 % <SyntaxCode>
304 out = cohere(at1,at5);
305 % </SyntaxCode>
306 stest = true;
307 catch err
308 disp(err.message)
309 stest = false;
310 end
311
312 % <AlgoDescription>
313 %
314 % 1) Check that the number of elements in 'out' is equal to 1.
315 % 2) Check that each output AO contains the correct data.
316 %
317 % </AlgoDescription>
318
319 atest = true;
320 if stest
321 % <AlgoCode>
322 % Check we have the correct number of outputs
323 if numel(out) ~= 1, atest = false; end
324
325 TOL = 1e-13;
326
327 % Get shortest vector
328 lmin = min([length(at1.y), length(at5.y)]);
329 % Set Nfft
330 Nfft = lmin;
331 % Get default window
332 if strcmpi(defaultWinType, 'kaiser')
333 win = specwin(defaultWinType, Nfft, find(ao.getInfo('cohere').plists, 'psll'));
334 else
335 win = specwin(defaultWinType, Nfft);
336 end
337 % Compute magnitude squared coherence estimate with MATLAB
338 % out: at1->at5
339 [cxy, f] = mscohere(at1.y(1:lmin), at5.y(1:lmin), win.win, Nfft/2, Nfft, at1.fs);
340 if any(abs(out.y-cxy > TOL)), atest = false; end
341 if any(abs(out.x-f > TOL)), atest = false; end
342 % </AlgoCode>
343 else
344 atest = false;
345 end
346
347 % Return a result structure
348 result = utp_prepare_result(atest, stest, dbstack, mfilename);
349 end % END UTP_04
350
351 %% UTP_05
352
353 % <TestDescription>
354 %
355 % Test that the cohere method doesn't work with an input of matrices
356 % and vectors and single AOs.
357 %
358 % </TestDescription>
359 function result = utp_05
360
361 % <SyntaxDescription>
362 %
363 % Test that the cohere method doesn't work with an input of matrices
364 % and vectors and single AOs.
365 %
366 % </SyntaxDescription>
367
368 try
369 % <SyntaxCode>
370 out = cohere([at5 at6], [at5 at1; at6 at1], at6);
371 stest = false;
372 % </SyntaxCode>
373 catch err
374 stest = true;
375 end
376
377 % <AlgoDescription>
378 %
379 % 1) Nothing to check
380 %
381 % </AlgoDescription>
382
383 atest = true;
384 if stest
385 % <AlgoCode>
386 % </AlgoCode>
387 else
388 atest = false;
389 end
390
391 % Return a result structure
392 result = utp_prepare_result(atest, stest, dbstack, mfilename);
393 end % END UTP_05
394
395 %% UTP_06
396
397 % <TestDescription>
398 %
399 % Tests that the cohere method properly applies history.
400 %
401 % </TestDescription>
402 function result = utp_06
403
404 % <SyntaxDescription>
405 %
406 % Test that the result of applying the cohere method can be processed back
407 % to an m-file.
408 %
409 % </SyntaxDescription>
410
411 try
412 % <SyntaxCode>
413 out = cohere(at5,at6);
414 mout = rebuild(out);
415 % </SyntaxCode>
416 stest = true;
417 catch err
418 disp(err.message)
419 stest = false;
420 end
421
422 % <AlgoDescription>
423 %
424 % 1) Check that the last entry in the history of 'out' corresponds to
425 % 'cohere'.
426 % 2) Check that the re-built object is the same as 'out'.
427 %
428 % </AlgoDescription>
429
430 atest = true;
431 if stest
432 % <AlgoCode>
433 % Check the last step in the history of 'out'
434 if ~strcmp(out.hist.methodInfo.mname, 'cohere'), atest = false; end
435 % Check the re-built object
436 if ~eq(mout, out, ple2), atest = false; end
437 % </AlgoCode>
438 else
439 atest = false;
440 end
441
442 % Return a result structure
443 result = utp_prepare_result(atest, stest, dbstack, mfilename);
444 end % END UTP_06
445
446 %% UTP_07
447
448 % <TestDescription>
449 %
450 % Tests that the cohere method can not modify the input AO.
451 %
452 % </TestDescription>
453 function result = utp_07
454
455 % <SyntaxDescription>
456 %
457 % Test that the cohere method can not modify the input AO.
458 % The method must throw an error for the modifier call.
459 %
460 % </SyntaxDescription>
461
462 try
463 % <SyntaxCode>
464 % copy at1 to work with
465 ain = ao(at1);
466 % modify ain
467 ain.cohere(at5);
468 % </SyntaxCode>
469 stest = false;
470 catch err
471 stest = true;
472 end
473
474 % <AlgoDescription>
475 %
476 % 1) Nothing to check.
477 %
478 % </AlgoDescription>
479
480 atest = true;
481 if stest
482 % <AlgoCode>
483 % </AlgoCode>
484 else
485 atest = false;
486 end
487
488 % Return a result structure
489 result = utp_prepare_result(atest, stest, dbstack, mfilename);
490 end % END UTP_07
491
492 %% UTP_08
493
494 % <TestDescription>
495 %
496 % Test the shape of the output.
497 %
498 % </TestDescription>
499 function result = utp_08
500
501 % <SyntaxDescription>
502 %
503 % Test that the cohere method keeps the data shape of the input object. The
504 % input AO must be an AO with row data and an AO with column data.
505 %
506 % </SyntaxDescription>
507
508 try
509 % <SyntaxCode>
510 out1 = cohere(at5, at6);
511 out2 = cohere(at6, at5);
512 % </SyntaxCode>
513 stest = true;
514 catch err
515 disp(err.message)
516 stest = false;
517 end
518
519 % <AlgoDescription>
520 %
521 % 1) Check that the shpe of the output data doesn't change.
522 %
523 % </AlgoDescription>
524
525 atest = true;
526 if stest
527 % <AlgoCode>
528 % Check the shape of the output data
529 if size(out1.data.y, 2) ~= 1, atest = false; end
530 if size(out2.data.y, 1) ~= 1, atest = false; end
531 % </AlgoCode>
532 else
533 atest = false;
534 end
535
536 % Return a result structure
537 result = utp_prepare_result(atest, stest, dbstack, mfilename);
538 end % END UTP_08
539
540 %% UTP_09
541
542 % <TestDescription>
543 %
544 % Check that the cohere method pass back the output objects to a list of
545 % output variables or to a single variable.
546 %
547 % </TestDescription>
548 function result = utp_09
549
550 % <SyntaxDescription>
551 %
552 % This test is not longer necessary because the cohere method pass back
553 % always only one object.
554 %
555 % </SyntaxDescription>
556
557 try
558 % <SyntaxCode>
559 % </SyntaxCode>
560 stest = true;
561 catch err
562 disp(err.message)
563 stest = false;
564 end
565
566 % <AlgoDescription>
567 %
568 % 1) Nothing to check.
569 %
570 % </AlgoDescription>
571
572 atest = true;
573 if stest
574 % <AlgoCode>
575 % </AlgoCode>
576 else
577 atest = false;
578 end
579
580 % Return a result structure
581 result = utp_prepare_result(atest, stest, dbstack, mfilename);
582 end % END UTP_09
583
584 %% UTP_10
585
586 % <TestDescription>
587 %
588 % Tests that the cohere method agrees with MATLAB's mscohere when
589 % configured to use the same parameters.
590 %
591 % </TestDescription>
592 function result = utp_10
593
594 % <SyntaxDescription>
595 %
596 % Test that applying cohere works on two AOs.
597 %
598 % </SyntaxDescription>
599
600 try
601 % <SyntaxCode>
602 % Construct two test AOs
603 nsecs = 10;
604 fs = 1000;
605 pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');
606 a1 = ao(pl); a2 = ao(pl);
607 % Filter one time-series
608 f2 = miir(plist('type', 'bandpass', 'fs', fs, 'order', 3, 'fc', [50 250]));
609 a1f = filter(a1, plist('filter', f2));
610 % make some cross-power
611 a4 = a1f+a2; a4.setName;
612 % Compute coherence
613 Nfft = 2*fs;
614 win = specwin('Hanning', Nfft);
615 pl = plist('Nfft', Nfft, 'Win', win.type, 'order', -1, 'type', 'MS');
616 out = cohere(a4,a1,pl);
617 % </SyntaxCode>
618 stest = true;
619 catch err
620 disp(err.message)
621 stest = false;
622 end
623
624 % <AlgoDescription>
625 %
626 % 1) Check that output agrees with the output of MATLAB's mscohere.
627 % 2) Check that the shape of the output data is equal to the input data
628 %
629 % </AlgoDescription>
630
631 atest = true;
632 if stest
633 % <AlgoCode>
634 % Compute coherence using MATLAB's cohere
635 [cxy, f] = mscohere(a4.y, a1.y, win.win, Nfft/2, Nfft, a1.fs);
636 if ne(cxy(:), out.y), atest = false; end
637 if ne(f, out.x), atest = false; end
638 if ne(out, out, ple2), atest = false; end
639 % Check the data shape
640 if size(a4.y,1) == 1
641 if size(out.y,1) ~= 1, atest = false; end
642 else
643 if size(out.y,2) ~= 1, atest = false; end
644 end
645 % </AlgoCode>
646 else
647 atest = false;
648 end
649
650 % Return a result structure
651 result = utp_prepare_result(atest, stest, dbstack, mfilename);
652 end % END UTP_10
653
654
655 %% UTP_12
656
657 % <TestDescription>
658 %
659 % Tests symmetry properties of complex-coherence:
660 % 1) white noise produced from normal pdf, with a given mean value and
661 % sigma (distribution's 1st and 2nd orders)
662 % 2) white noise produced from normal pdf, with a given mean value and
663 % sigma (distribution's 1st and 2nd orders)
664 % 3) complex coherence of the white noise series
665 % 4) compare C(x,y) with conj(C(y,x))
666 % 5) compare C(x,x) and C(y,y) with 1
667 %
668
669 % </TestDescription>
670 function result = utp_12
671
672 % <SyntaxDescription>
673 %
674 % 1) Prepare the test tsdata:
675 % white noise from normal distribution + offset
676 % 2) Assign a random unit
677 % 3) Prepare the test tsdata:
678 % white noise from normal distribution + offset
679 % 4) Assign a random unit
680 % 5) complex coherence of the white noise
681 %
682 % </SyntaxDescription>
683
684 % <SyntaxCode>
685 try
686
687 % Array of parameters to pick from
688 fs_list = [0.1;1;10];
689 nsecs_list = [100:100:10000]';
690 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
691 mu_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
692
693 % Build time-series test data
694
695 % Picks the values at random from the list
696 fs = utils.math.randelement(fs_list, 1);
697 nsecs = utils.math.randelement(nsecs_list, 1);
698 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
699 mu_distr = utils.math.randelement(mu_distr_list, 1);
700 f = [1:5] / 100 * fs;
701 A = sigma_distr + sigma_distr*rand(1,1);
702 phi = 0 + 2*pi*rand(1,1);
703
704 % White noise
705 type = 'Normal';
706 a_n1 = ao(plist('waveform', 'noise', ...
707 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
708 a_n2 = ao(plist('waveform', 'noise', ...
709 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
710 a_const = ao(mu_distr);
711 a_wave = ao(plist('waveform', 'sine-wave', ...
712 'fs', fs, 'nsecs', nsecs, 'f', f, 'A', A, 'phi', phi));
713 a_1 = a_n1 + a_const + a_wave;
714 a_2 = a_n2 + a_wave;
715
716 % Set units and prefix from those supported
717 unit_list = unit.supportedUnits;
718 % remove the first empty unit '' from the list, because then is it
719 % possible that we add a prefix to an empty unit
720 unit_list = unit_list(2:end);
721 prefix_list = unit.supportedPrefixes;
722 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
723 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
724
725 % Evaluate the complex coherence of the time-series data
726 win_list = specwin.getTypes;
727 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
728 win_type = win_type{1};
729 if strcmp(win_type, 'Kaiser')
730 win = specwin(win_type, 1, find(ao.getInfo('psd').plists, 'psll'));
731 else
732 win = specwin(win_type, 1);
733 end
734 olap = win.rov;
735 detrend = 0;
736 scale_type = 'C';
737 n_pts = nsecs*fs/10;
738
739 C12 = cohere(a_1, a_2, ...
740 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
741 C21 = cohere(a_2, a_1, ...
742 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
743 C21_cc = conj(C21);
744 C11 = cohere(a_1, a_1, ...
745 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
746 C22 = cohere(a_2, a_2, ...
747 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
748 stest = true;
749
750 catch err
751 disp(err.message)
752 stest = false;
753 end
754 % </SyntaxCode>
755
756 % <AlgoDescription>
757 %
758 % 1) Check that C(x,y) equals conj(C(y,x))
759 % 2) Check that C(x,x) equals 1
760 % 2) Check that C(y,y) equals 1
761
762 % </AlgoDescription>
763
764 % <AlgoCode>
765 atest = true;
766
767 if stest
768 if ~eq(C12.data, C21_cc.data, 'dy') || ...
769 ~isequal(C11.y, ones(size(C11.y))) || ...
770 ~isequal(C22.y, ones(size(C22.y)))
771 atest = false;
772 end
773 else
774 atest = false;
775 end
776 % </AlgoCode>
777
778 % Return a result structure
779 result = utp_prepare_result(atest, stest, dbstack, mfilename);
780 end % END UTP_12
781
782 %% UTP_13
783
784 % <TestDescription>
785 %
786 % Tests symmetry properties of complex-coherence:
787 % 1) white noise produced from normal pdf, with a given mean value and
788 % sigma (distribution's 1st and 2nd orders)
789 % 2) white noise produced from normal pdf, with a given mean value and
790 % sigma (distribution's 1st and 2nd orders)
791 % 3) magnitude-squared coherence of the white noise series
792 % 4) compare C(x,y) with C(y,x)
793 % 5) compare C(x,x) and C(y,y) with 1
794 %
795
796 % </TestDescription>
797 function result = utp_13
798
799 % <SyntaxDescription>
800 %
801 % 1) Prepare the test tsdata:
802 % white noise from normal distribution + offset
803 % 2) Assign a random unit
804 % 3) Prepare the test tsdata:
805 % white noise from normal distribution + offset
806 % 4) Assign a random unit
807 % 5) magnitude-squared coherence of the white noise
808 %
809 % </SyntaxDescription>
810
811 % <SyntaxCode>
812 try
813
814 % Array of parameters to pick from
815 fs_list = [0.1;1;10];
816 nsecs_list = [100:100:10000]';
817 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
818 mu_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
819
820 % Build time-series test data
821
822 % Picks the values at random from the list
823 fs = utils.math.randelement(fs_list, 1);
824 nsecs = utils.math.randelement(nsecs_list, 1);
825 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
826 mu_distr = utils.math.randelement(mu_distr_list, 1);
827 f = [1:5] / 100 * fs;
828 A = sigma_distr + sigma_distr*rand(1,1);
829 phi = 0 + 2*pi*rand(1,1);
830
831 % White noise
832 type = 'Normal';
833 a_n1 = ao(plist('waveform', 'noise', ...
834 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
835 a_n2 = ao(plist('waveform', 'noise', ...
836 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
837 a_const = ao(mu_distr);
838 a_wave = ao(plist('waveform', 'sine-wave', ...
839 'fs', fs, 'nsecs', nsecs, 'f', f, 'A', A, 'phi', phi));
840 a_1 = a_n1 + a_const + a_wave;
841 a_2 = a_n2 + a_wave;
842
843 % Set units and prefix from those supported
844 unit_list = unit.supportedUnits;
845 % remove the first empty unit '' from the list, because then is it
846 % possible that we add a prefix to an empty unit
847 unit_list = unit_list(2:end);
848 prefix_list = unit.supportedPrefixes;
849 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
850 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
851
852 % Evaluate the magnitude-squared coherence of the time-series data
853 win_list = specwin.getTypes;
854 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
855 win_type = win_type{1};
856 if strcmp(win_type, 'Kaiser')
857 win = specwin(win_type, 1, find(ao.getInfo('psd').plists, 'psll'));
858 else
859 win = specwin(win_type, 1);
860 end
861 olap = win.rov;
862 detrend = 0;
863 scale_type = 'MS';
864 n_pts = nsecs*fs/10;
865
866 C12 = cohere(a_1, a_2, ...
867 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
868 C21 = cohere(a_2, a_1, ...
869 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
870 C11 = cohere(a_1, a_1, ...
871 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
872 C22 = cohere(a_2, a_2, ...
873 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
874 stest = true;
875
876 catch err
877 disp(err.message)
878 stest = false;
879 end
880 % </SyntaxCode>
881
882 % <AlgoDescription>
883 %
884 % 1) Check that C(x,y) equals C(y,x)
885 % 1) Check that C(x,x) equals 1
886 % 1) Check that C(y,y) equals 1
887
888 % </AlgoDescription>
889
890 % <AlgoCode>
891 atest = true;
892
893 if stest
894 if ~isequal(C12.data, C21.data) || ...
895 ~isequal(C11.y, ones(size(C11.y))) ...
896 || ~isequal(C22.y, ones(size(C22.y)))
897 atest = false;
898 end
899 if atest == false
900 fs
901 nsecs
902 sigma_distr
903 mu_distr
904 f
905 A
906 phi
907 end
908 else
909 atest = false;
910 end
911 % </AlgoCode>
912
913 % Return a result structure
914 result = utp_prepare_result(atest, stest, dbstack, mfilename);
915 end % END UTP_13
916
917 %% UTP_14
918
919 % <TestDescription>
920 %
921 % Tests symmetry properties of complex-coherence:
922 % 1) white noise produced from normal pdf, with a given mean value and
923 % sigma (distribution's 1st and 2nd orders)
924 % 2) white noise produced from normal pdf, with a given mean value and
925 % sigma (distribution's 1st and 2nd orders)
926 % 3) complex coherence of the combination of white noise series
927 % 4) compare C(x,y) with 1
928 %
929
930 % </TestDescription>
931 function result = utp_14
932
933 % <SyntaxDescription>
934 %
935 % 1) Prepare the test tsdata:
936 % white noise from normal distribution + offset
937 % 2) Assign a random unit
938 % 3) Prepare the test tsdata:
939 % white noise from normal distribution + offset
940 % 4) Assign a random unit
941 % 5) complex coherence of the combination of noise
942 %
943 % </SyntaxDescription>
944
945 % <SyntaxCode>
946 try
947
948 % Array of parameters to pick from
949 fs_list = [0.1;1;10];
950 nsecs_list = [100:100:10000]';
951 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
952 mu_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
953
954 % Build time-series test data
955
956 % Picks the values at random from the list
957 fs = utils.math.randelement(fs_list, 1);
958 nsecs = utils.math.randelement(nsecs_list, 1);
959 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
960 mu_distr = utils.math.randelement(mu_distr_list, 1);
961 f = [1:5] / 100 * fs;
962 A = sigma_distr + sigma_distr*rand(1,1);
963 phi = 0 + 2*pi*rand(1,1);
964
965 % White noise
966 type = 'Normal';
967 a_n = ao(plist('waveform', 'noise', ...
968 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
969 a_const = ao(mu_distr);
970 % Sinusoidal signal
971 a_wave = ao(plist('waveform', 'sine-wave', ...
972 'fs', fs, 'nsecs', nsecs, 'f', f, 'A', A, 'phi', phi));
973 a_1 = a_n + a_wave;
974 % Linear combination (totally correlated time series)
975 a_2 = a_1 + a_const;
976
977 % Set units and prefix from those supported
978 unit_list = unit.supportedUnits;
979 % remove the first empty unit '' from the list, because then is it
980 % possible that we add a prefix to an empty unit
981 unit_list = unit_list(2:end);
982 prefix_list = unit.supportedPrefixes;
983 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
984 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
985
986 % Evaluate the complex coherence of the time-series data
987 win_list = specwin.getTypes;
988 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
989 win_type = win_type{1};
990 if strcmp(win_type, 'Kaiser')
991 win = specwin(win_type, 1, find(ao.getInfo('psd').plists, 'psll'));
992 else
993 win = specwin(win_type, 1);
994 end
995 olap = win.rov;
996 detrend = 0;
997 scale_type = 'C';
998 n_pts = nsecs*fs/10;
999
1000 C = cohere(a_1, a_2, ...
1001 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
1002 stest = true;
1003
1004 catch err
1005 disp(err.message)
1006 stest = false;
1007 end
1008 % </SyntaxCode>
1009
1010 % <AlgoDescription>
1011 %
1012 % 1) Check that the complex coherence equals 1
1013
1014 % </AlgoDescription>
1015
1016 % <AlgoCode>
1017 atest = true;
1018 TOL = 1e-12;
1019
1020 if stest
1021 if any(abs((C.y - 1)) > TOL)
1022 atest = false;
1023 end
1024 else
1025 atest = false;
1026 end
1027 % </AlgoCode>
1028
1029 % Return a result structure
1030 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1031 end % END UTP_14
1032
1033 %% UTP_15
1034
1035 % <TestDescription>
1036 %
1037 % Tests symmetry properties of complex-coherence:
1038 % 1) white noise produced from normal pdf, with a given mean value and
1039 % sigma (distribution's 1st and 2nd orders)
1040 % 2) white noise produced from normal pdf, with a given mean value and
1041 % sigma (distribution's 1st and 2nd orders)
1042 % 3) magnitude-squared coherence of the combination of white noise series
1043 % 4) compare C(x,y) with 1
1044 %
1045
1046 % </TestDescription>
1047 function result = utp_15
1048
1049 % <SyntaxDescription>
1050 %
1051 % 1) Prepare the test tsdata:
1052 % white noise from normal distribution + offset
1053 % 2) Assign a random unit
1054 % 3) Prepare the test tsdata:
1055 % white noise from normal distribution + offset
1056 % 4) Assign a random unit
1057 % 5) magnitude-squared coherence of the combination of noise
1058 %
1059 % </SyntaxDescription>
1060
1061 % <SyntaxCode>
1062 try
1063
1064 % Array of parameters to pick from
1065 fs_list = [0.1;1;10];
1066 nsecs_list = [100:100:10000]';
1067 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1068 mu_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1069
1070 % Build time-series test data
1071
1072 % Picks the values at random from the list
1073 fs = utils.math.randelement(fs_list, 1);
1074 nsecs = utils.math.randelement(nsecs_list, 1);
1075 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
1076 mu_distr = utils.math.randelement(mu_distr_list, 1);
1077 f = [1:5] / 100 * fs;
1078 A = sigma_distr + sigma_distr*rand(1,1);
1079 phi = 0 + 2*pi*rand(1,1);
1080
1081 % White noise
1082 type = 'Normal';
1083 a_n = ao(plist('waveform', 'noise', ...
1084 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1085 a_const = ao(mu_distr);
1086 % Sinusoidal signal
1087 a_wave = ao(plist('waveform', 'sine-wave', ...
1088 'fs', fs, 'nsecs', nsecs, 'f', f, 'A', A, 'phi', phi));
1089 a_1 = a_n + a_wave;
1090 % Linear combination (totally correlated time series)
1091 a_2 = a_1 + a_const;
1092
1093 % Set units and prefix from those supported
1094 unit_list = unit.supportedUnits;
1095 % remove the first empty unit '' from the list, because then is it
1096 % possible that we add a prefix to an empty unit
1097 unit_list = unit_list(2:end);
1098 prefix_list = unit.supportedPrefixes;
1099 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1100 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1101
1102 % Evaluate the complex coherence of the time-series data
1103 win_list = specwin.getTypes;
1104 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
1105 win_type = win_type{1};
1106 if strcmp(win_type, 'Kaiser')
1107 win = specwin(win_type, 1, find(ao.getInfo('psd').plists, 'psll'));
1108 else
1109 win = specwin(win_type, 1);
1110 end
1111 olap = win.rov;
1112 detrend = 0;
1113 scale_type = 'MS';
1114 n_pts = nsecs*fs/10;
1115
1116 C = cohere(a_1, a_2, ...
1117 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
1118 stest = true;
1119
1120 catch err
1121 disp(err.message)
1122 stest = false;
1123 end
1124 % </SyntaxCode>
1125
1126 % <AlgoDescription>
1127 %
1128 % 1) Check that the magnitude-squared coherence equals 1
1129
1130 % </AlgoDescription>
1131
1132 % <AlgoCode>
1133 atest = true;
1134
1135 if stest
1136 if ~eq(C.y, ones(size(C.y)))
1137 atest = false;
1138 end
1139 else
1140 atest = false;
1141 end
1142 % </AlgoCode>
1143
1144 % Return a result structure
1145 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1146 end % END UTP_15
1147
1148 %% UTP_16
1149
1150 % <TestDescription>
1151 %
1152 % Tests symmetry properties of complex-coherence:
1153 % 1) white noise produced from normal pdf, with a given mean value and
1154 % sigma (distribution's 1st and 2nd orders)
1155 % 2) white noise produced from normal pdf, with a given mean value and
1156 % sigma (distribution's 1st and 2nd orders)
1157 % 3) magnitude-squared coherence M of the combination of white noise series
1158 % 4) complex coherence C of the combination of white noise series
1159 % 5) compare abs(C)^2 with M
1160 %
1161
1162 % </TestDescription>
1163 function result = utp_16
1164
1165 % <SyntaxDescription>
1166 %
1167 % 1) Prepare the test tsdata:
1168 % white noise from normal distribution + offset
1169 % 2) Assign a random unit
1170 % 3) Prepare the test tsdata:
1171 % white noise from normal distribution + offset
1172 % 4) Assign a random unit
1173 % 5) magnitude-squared coherence of the combination of noise
1174 % 6) complex coherence of the combination of noise
1175 %
1176 % </SyntaxDescription>
1177
1178 % <SyntaxCode>
1179 try
1180
1181 % Array of parameters to pick from
1182 fs_list = [0.1;1;10];
1183 nsecs_list = [100:100:10000]';
1184 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1185 mu_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1186
1187 % Build time-series test data
1188
1189 % Picks the values at random from the list
1190 fs = utils.math.randelement(fs_list, 1);
1191 nsecs = utils.math.randelement(nsecs_list, 1);
1192 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
1193 mu_distr = utils.math.randelement(mu_distr_list, 1);
1194 f = [1:5] / 100 * fs;
1195 A = sigma_distr + sigma_distr*rand(1,1);
1196 phi = 0 + 2*pi*rand(1,1);
1197
1198 % White noise
1199 type = 'Normal';
1200 a_n = ao(plist('waveform', 'noise', ...
1201 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1202 a_const = ao(mu_distr);
1203 % Sinusoidal signal
1204 a_wave = ao(plist('waveform', 'sine-wave', ...
1205 'fs', fs, 'nsecs', nsecs, 'f', f, 'A', A, 'phi', phi));
1206 a_1 = a_n + a_wave;
1207 % Linear combination (totally correlated time series)
1208 a_2 = a_1 + a_const;
1209
1210 % Set units and prefix from those supported
1211 unit_list = unit.supportedUnits;
1212 % remove the first empty unit '' from the list, because then is it
1213 % possible that we add a prefix to an empty unit
1214 unit_list = unit_list(2:end);
1215 prefix_list = unit.supportedPrefixes;
1216 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1217 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1218
1219 % Evaluate the complex coherence of the time-series data
1220 win_list = specwin.getTypes;
1221 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
1222 win_type = win_type{1};
1223 if strcmp(win_type, 'Kaiser')
1224 win = specwin(win_type, 1, find(ao.getInfo('psd').plists, 'psll'));
1225 else
1226 win = specwin(win_type, 1);
1227 end
1228 olap = win.rov;
1229 detrend = 0;
1230 n_pts = nsecs*fs/10;
1231
1232 M = cohere(a_1, a_2, ...
1233 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', 'MS'));
1234 C = cohere(a_1, a_2, ...
1235 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', 'C'));
1236 stest = true;
1237
1238 catch err
1239 disp(err.message)
1240 stest = false;
1241 end
1242 % </SyntaxCode>
1243
1244 % <AlgoDescription>
1245 %
1246 % 1) Check that the magnitude-squared coherence equals the square
1247 % modulus of the complex coherence
1248
1249 % </AlgoDescription>
1250
1251 % <AlgoCode>
1252 atest = true;
1253 TOL = 1e-15;
1254
1255 if stest
1256 if any(abs(M.y - abs(C.y).^2) > TOL)
1257 atest = false;
1258 end
1259 else
1260 atest = false;
1261 end
1262 % </AlgoCode>
1263
1264 % Return a result structure
1265 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1266 end % END UTP_16
1267
1268 %% UTP_17
1269
1270 % <TestDescription>
1271 %
1272 % Tests handling of units:
1273 % 1) white noise produced from normal pdf, with a given mean value and
1274 % sigma (distribution's 1st and 2nd orders)
1275 % 2) white noise produced from normal pdf, with a given mean value and
1276 % sigma (distribution's 1st and 2nd orders)
1277 % 3) complex coherence of the white noise series
1278 % 4) compares the units of the input and output
1279 %
1280
1281 % </TestDescription>
1282 function result = utp_17
1283
1284 % <SyntaxDescription>
1285 %
1286 % 1) Prepare the test tsdata:
1287 % white noise from normal distribution + offset
1288 % 2) Assign a random unit
1289 % 3) Prepare the test tsdata:
1290 % white noise from normal distribution + offset
1291 % 4) Assign a random unit
1292 % 5) complex cohere of the white noise
1293 %
1294 % </SyntaxDescription>
1295
1296 % <SyntaxCode>
1297 try
1298
1299 % Build time-series test data
1300 fs = 1;
1301 nsecs = 86400;
1302 sigma_distr_1 = 4.69e-12;
1303 mu_distr_1 = -5.11e-14;
1304 sigma_distr_2 = 6.04e-9;
1305 mu_distr_2 = 1.5e-10;
1306
1307 % White noise
1308 type = 'Normal';
1309
1310 a_n = ao(plist('waveform', 'noise', ...
1311 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1));
1312 a_const = ao(mu_distr_1);
1313 a_1 = a_n + a_const;
1314
1315 a_n = ao(plist('waveform', 'noise', ...
1316 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_2));
1317 a_const = ao(mu_distr_2);
1318 a_2 = a_n + a_const;
1319
1320 % Set units and prefix from those supported
1321 unit_list = unit.supportedUnits;
1322 % remove the first empty unit '' from the list, because then is it
1323 % possible that we add a prefix to an empty unit
1324 unit_list = unit_list(2:end);
1325 prefix_list = unit.supportedPrefixes;
1326 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1327 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1328
1329 % Evaluate the coherence of the time-series data
1330 win = specwin('BH92');
1331 olap = win.rov;
1332 detrend = 0;
1333 scale_type = 'C';
1334 n_pts = nsecs*fs/10;
1335
1336 C = cohere(a_1, a_2, ...
1337 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
1338
1339 stest = true;
1340
1341 catch err
1342 disp(err.message)
1343 stest = false;
1344 end
1345 % </SyntaxCode>
1346
1347 % <AlgoDescription>
1348 %
1349 % 1) Check that (complex coherence yunits) equals [1]
1350 % 2) Check that (complex coherence xunits) equals [Hz]
1351
1352 % </AlgoDescription>
1353
1354 % <AlgoCode>
1355 atest = true;
1356
1357 if stest
1358 if ~eq(C.yunits, unit('')) || ~eq(C.xunits, unit('Hz'))
1359 atest = false;
1360 end
1361 else
1362 atest = false;
1363 end
1364 % </AlgoCode>
1365
1366 % Return a result structure
1367 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1368 end % END UTP_17
1369
1370 %% UTP_18
1371
1372 % <TestDescription>
1373 %
1374 % Tests handling of units:
1375 % 1) white noise produced from normal pdf, with a given mean value and
1376 % sigma (distribution's 1st and 2nd orders)
1377 % 2) white noise produced from normal pdf, with a given mean value and
1378 % sigma (distribution's 1st and 2nd orders)
1379 % 3) magnitude-squared coherence of the white noise series
1380 % 4) compares the units of the input and output
1381 %
1382
1383 % </TestDescription>
1384 function result = utp_18
1385
1386 % <SyntaxDescription>
1387 %
1388 % 1) Prepare the test tsdata:
1389 % white noise from normal distribution + offset
1390 % 2) Assign a random unit
1391 % 3) Prepare the test tsdata:
1392 % white noise from normal distribution + offset
1393 % 4) Assign a random unit
1394 % 5) magnitude-squared cohere of the white noise
1395 %
1396 % </SyntaxDescription>
1397
1398 % <SyntaxCode>
1399 try
1400
1401 % Build time-series test data
1402 fs = 1;
1403 nsecs = 86400;
1404 sigma_distr_1 = 4.69e-12;
1405 mu_distr_1 = -5.11e-14;
1406 sigma_distr_2 = 6.04e-9;
1407 mu_distr_2 = 1.5e-10;
1408
1409 % White noise
1410 type = 'Normal';
1411
1412 a_n = ao(plist('waveform', 'noise', ...
1413 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1));
1414 a_const = ao(mu_distr_1);
1415 a_1 = a_n + a_const;
1416
1417 a_n = ao(plist('waveform', 'noise', ...
1418 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_2));
1419 a_const = ao(mu_distr_2);
1420 a_2 = a_n + a_const;
1421
1422 % Set units and prefix from those supported
1423 unit_list = unit.supportedUnits;
1424 % remove the first empty unit '' from the list, because then is it
1425 % possible that we add a prefix to an empty unit
1426 unit_list = unit_list(2:end);
1427 prefix_list = unit.supportedPrefixes;
1428 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1429 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1430
1431 % Evaluate the coherence of the time-series data
1432 win = specwin('BH92');
1433 olap = win.rov;
1434 detrend = 0;
1435 scale_type = 'MS';
1436 n_pts = nsecs*fs/10;
1437
1438 C = cohere(a_1, a_2, ...
1439 plist('Win', win.type, 'olap', olap, 'Nfft', n_pts, 'order', detrend,'type', scale_type));
1440
1441 stest = true;
1442
1443 catch err
1444 disp(err.message)
1445 stest = false;
1446 end
1447 % </SyntaxCode>
1448
1449 % <AlgoDescription>
1450 %
1451 % 1) Check that (magnitude-squared coherence yunits) equals [1]
1452 % 2) Check that (magnitude-squared coherence xunits) equals [Hz]
1453
1454 % </AlgoDescription>
1455
1456 % <AlgoCode>
1457 atest = true;
1458
1459 if stest
1460 if ~eq(C.yunits, unit('')) || ~eq(C.xunits, unit('Hz'))
1461 atest = false;
1462 end
1463 else
1464 atest = false;
1465 end
1466 % </AlgoCode>
1467
1468 % Return a result structure
1469 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1470 end % END UTP_18
1471
1472 %% UTP_19
1473
1474 % <TestDescription>
1475 %
1476 % Tests that differently sized data sets are treated properly
1477 %
1478 % </TestDescription>
1479 function result = utp_19
1480
1481 % <SyntaxDescription>
1482 %
1483 % Test that applying cohere works on two AOs.
1484 %
1485 % </SyntaxDescription>
1486
1487 try
1488 % <SyntaxCode>
1489 % Construct two test AOs
1490 nsecs = [10000:1:20000];
1491 fs = 1;
1492 pl = plist('fs', fs, 'tsfcn', 'randn(size(t))');
1493 a1 = ao(pl.pset('nsecs', utils.math.randelement(nsecs, 1)));
1494 a2 = ao(pl.pset('nsecs', utils.math.randelement(nsecs, 1)));
1495 len_1 = a1.len;
1496 len_2 = a2.len;
1497 % Filter one time-series
1498 f2 = miir(plist('type', 'bandpass', 'fs', fs, 'order', 3, 'fc', [.050 .25]));
1499 a1f = filter(a1, plist('filter', f2));
1500 % Compute cohere
1501 Nfft = -1;
1502 win = 'Hanning';
1503 pl = plist('Nfft', Nfft, 'Win', win, 'order', -1);
1504 out = cohere(a2,a1f,pl);
1505 % </SyntaxCode>
1506 stest = true;
1507 catch err
1508 disp(err.message)
1509 stest = false;
1510 end
1511
1512 % <AlgoDescription>
1513 %
1514 % 1) Check that cohere used the length of the shortest ao.
1515 %
1516 % </AlgoDescription>
1517
1518 atest = true;
1519 if stest
1520 % <AlgoCode>
1521 % Compare the nfft with the length of the input data
1522
1523 if out.x(2) ~= 1/min(len_1,len_2)
1524 atest = false;
1525 end
1526 % </AlgoCode>
1527 else
1528 atest = false;
1529 end
1530
1531 % Return a result structure
1532 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1533 end % END UTP_19
1534
1535 %% UTP_20
1536
1537 % <TestDescription>
1538 %
1539 % Tests that applying a single window the coherence is 1
1540 %
1541 % </TestDescription>
1542 function result = utp_20
1543
1544 % <SyntaxDescription>
1545 %
1546 % Test that applying cohere works on two AOs.
1547 %
1548 % </SyntaxDescription>
1549
1550 try
1551 % <SyntaxCode>
1552 % Construct two test AOs
1553 nsecs = [10000:100:20000];
1554 fs = 1;
1555 pl = plist('fs', fs, 'tsfcn', 'randn(size(t))');
1556 a1 = ao(pl.pset('nsecs', utils.math.randelement(nsecs, 1)));
1557 a2 = ao(pl.pset('nsecs', utils.math.randelement(nsecs, 1)));
1558 % Filter one time-series
1559 f2 = miir(plist('type', 'bandpass', 'fs', fs, 'order', 3, 'fc', [.050 .25]));
1560 a1f = filter(a1, plist('filter', f2));
1561 % Compute cohere
1562 Nfft = -1;
1563 win = 'Hanning';
1564 pl = plist('Nfft', Nfft, 'Win', win, 'order', -1);
1565 out_c = cohere(a2, a1f, pl.pset('type', 'C'));
1566 out_ms = cohere(a2, a1f, pl.pset('type', 'MS'));
1567 % </SyntaxCode>
1568 stest = true;
1569 catch err
1570 disp(err.message)
1571 stest = false;
1572 end
1573
1574 % <AlgoDescription>
1575 %
1576 % 1) Check that the calculated cohere is 1
1577 %
1578 % </AlgoDescription>
1579
1580 atest = true;
1581 TOL = 1e-12;
1582 if stest
1583 % <AlgoCode>
1584 % Compare the calculated cohere with 1
1585
1586 if any(abs(abs(out_c.y) - 1) > TOL)
1587 atest = false;
1588 end
1589 if any(abs(abs(out_ms.y) - 1) > TOL)
1590 atest = false;
1591 end
1592 % </AlgoCode>
1593 else
1594 atest = false;
1595 end
1596
1597 % Return a result structure
1598 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1599 end % END UTP_20
1600
1601 %% UTP_21
1602
1603 % <TestDescription>
1604 %
1605 % Tests the possibility to set the number of averages rather than setting the Nfft:
1606 % 1) white noise produced from normal pdf, with:
1607 % a given mean value and sigma (distribution's 1st and 2nd order)
1608 % 2) cohere of the noise, without detrending, random window, set number of
1609 % averages
1610 % 3) check the effective number of averages
1611 %
1612
1613 % </TestDescription>
1614 function result = utp_21
1615
1616 % <SyntaxDescription>
1617 %
1618 % 1) Prepare the test tsdata:
1619 % white noise from normal distribution + offset
1620 % 2) cohere of the noise, without detrending, random window, set number of
1621 % averages
1622 %
1623 % </SyntaxDescription>
1624
1625 % <SyntaxCode>
1626 try
1627 % Array of parameters to pick from
1628 fs_list = [0.1;1;2;5;10];
1629 nsecs_list = [2000:1000:10000]';
1630 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1631 trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]';
1632
1633 % Build time-series test data
1634
1635 % Picks the values at random from the list
1636 fs = utils.math.randelement(fs_list, 1);
1637 nsecs = utils.math.randelement(nsecs_list, 1);
1638 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
1639 trend_0 = utils.math.randelement(trend_0_list, 1);
1640
1641 % White noise
1642 type = 'Normal';
1643 a_n1 = ao(plist('waveform', 'noise', ...
1644 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1645 a_n2 = ao(plist('waveform', 'noise', ...
1646 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1647
1648 % Constant signal
1649 a_c = ao(trend_0);
1650
1651 % Total signals
1652 a1 = a_n1 + a_c;
1653 a2 = a_n2 + a_c;
1654
1655 % Evaluate the complex coherence of the white noise time-series data
1656 win_list = specwin.getTypes;
1657 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
1658 win_type = win_type{1};
1659 switch win_type
1660 case 'Kaiser'
1661 win = specwin(win_type, 1, find(ao.getInfo('psd').plists, 'psll'));
1662 otherwise
1663 win = specwin(win_type, 1);
1664 end
1665
1666 olap = win.rov;
1667 detrend = 0;
1668 n_pts = -1;
1669 scale_type = 'C';
1670 navs = utils.math.randelement([1:100],1);
1671
1672 % Evaluates the coherence asking for the number of averages
1673 C = cohere(a1, a2, plist('Win', win.type, 'olap', olap, ...
1674 'Nfft', n_pts, 'order', detrend, 'type', scale_type, 'navs', navs));
1675
1676 stest = true;
1677
1678 catch err
1679 disp(err.message)
1680 stest = false;
1681 end
1682 % </SyntaxCode>
1683
1684 % <AlgoDescription>
1685 %
1686 % 1) Check that calculated navs are identical to those requested
1687 %
1688 % </AlgoDescription>
1689
1690 % <AlgoCode>
1691 atest = true;
1692
1693 if stest
1694 % Compare the navs written in the output object with the requested one
1695 if ne(navs, C.data.navs)
1696 if ne(find(C.hist.plistUsed, 'navs'), C.data.navs)
1697 atest = false;
1698 end
1699 end
1700 else
1701 atest = false;
1702 end
1703 % </AlgoCode>
1704
1705 % Return a result structure
1706 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1707 end % END UTP_21
1708
1709 %% UTP_22
1710
1711 % <TestDescription>
1712 %
1713 % Tests the possibility to set the number of averages rather than setting the Nfft:
1714 % 1) white noise produced from uniform pdf, with:
1715 % a given mean value and sigma (distribution's 1st and 2nd order)
1716 % 2) cohere of the noise, without detrending, random window, random navs
1717 % 3) get the number of averages
1718 % 4) get the nfft used
1719 % 5) run cohere again, with the nfft used
1720 % 6) compare the calculated objects
1721 %
1722
1723 % </TestDescription>
1724 function result = utp_22
1725
1726 % <SyntaxDescription>
1727 %
1728 % 1) white noise produced from uniform pdf, with:
1729 % a given mean value and sigma (distribution's 1st and 2nd order)
1730 % 2) cohere of the noise, without detrending, random window, random navs
1731 % 3) get the number of averages
1732 % 4) get the nfft used
1733 % 5) run cohere again, with the nfft used
1734 %
1735 % </SyntaxDescription>
1736
1737 % <SyntaxCode>
1738 try
1739 % Array of parameters to pick from
1740 fs_list = [0.1;1;2;5;10];
1741 nsecs_list = [20 100 1000:1000:10000]';
1742 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1743 trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]';
1744
1745 % Build time-series test data
1746
1747 % Picks the values at random from the list
1748 fs = utils.math.randelement(fs_list, 1);
1749 nsecs = utils.math.randelement(nsecs_list, 1);
1750 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
1751 trend_0 = utils.math.randelement(trend_0_list, 1);
1752
1753 % White noise
1754 type = 'Uniform';
1755 a_n1 = ao(plist('waveform', 'noise', ...
1756 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1757 a_n2 = ao(plist('waveform', 'noise', ...
1758 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1759
1760 % Constant signal
1761 a_c = ao(trend_0);
1762
1763 % Total signals
1764 a1 = a_n1 + a_c;
1765 a2 = a_n2 + a_c;
1766
1767 % Evaluate the complex coherence of the white noise time-series data
1768 win_list = specwin.getTypes;
1769 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
1770 win_type = win_type{1};
1771 switch win_type
1772 case 'Kaiser'
1773 win = specwin(win_type, 1, find(ao.getInfo('psd').plists, 'psll'));
1774 otherwise
1775 win = specwin(win_type, 1);
1776 end
1777
1778 olap = win.rov;
1779 detrend = 0;
1780 scale_type = 'MS';
1781 navs = fix(utils.math.randelement(logspace(0,log10(max(0,a1.len/10)),50),1));
1782
1783 % Calculates the coherence asking for the number of averages
1784 C1 = cohere(a1, a2, plist('Win', win.type, 'olap', olap, ...
1785 'Nfft', -1, 'order', detrend, 'type', scale_type, ...
1786 'navs', navs));
1787
1788 % Calculates the coherence asking for the number of points just evaluated
1789 C2 = cohere(a1, a2, plist('Win', win.type, 'olap', olap, ...
1790 'Nfft', find(C1.hist.plistUsed, 'Nfft'), 'order', detrend, 'type', scale_type));
1791 stest = true;
1792
1793 catch err
1794 disp(err.message)
1795 stest = false;
1796 end
1797 % </SyntaxCode>
1798
1799 % <AlgoDescription>
1800 %
1801 % 1) Check that calculated objects C1 and C2 are identical
1802 %
1803 % </AlgoDescription>
1804
1805 % <AlgoCode>
1806 atest = true;
1807
1808 if stest
1809 % Compare the output objects
1810 if ne(C1,C2,ple3)
1811 atest = false;
1812 end
1813 else
1814 atest = false;
1815 end
1816 % </AlgoCode>
1817
1818 % Return a result structure
1819 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1820 end % END UTP_22
1821
1822 %% UTP_23
1823
1824 % <TestDescription>
1825 %
1826 % Tests the possibility to set the number of averages rather than setting the Nfft:
1827 % 1) white noise produced from normal pdf, with:
1828 % a given mean value and sigma (distribution's 1st and 2nd order)
1829 % 2) cohere of the noise, without detrending, random window, random navs
1830 % 3) get the number of averages
1831 % 4) get the nfft used
1832 % 5) run cohere again, with the nfft used
1833 % 6) compare navs, nfft, coheres
1834 %
1835
1836 % </TestDescription>
1837 function result = utp_23
1838
1839 % <SyntaxDescription>
1840 %
1841 % 1) white noise produced from normal pdf, with:
1842 % a given mean value and sigma (distribution's 1st and 2nd order)
1843 % 2) cohere of the noise, without detrending, random window, random navs
1844 % 3) get the number of averages
1845 % 4) get the nfft used
1846 % 5) run cohere again, with the nfft used
1847 % 6) run cohere again, with conflicting parameters, and verify it uses
1848 % nfft rather than navs
1849 %
1850 % </SyntaxDescription>
1851
1852 % <SyntaxCode>
1853 try
1854 % Array of parameters to pick from
1855 fs_list = [0.1;1;2;5;10];
1856 nsecs_list = [1000:1000:10000]';
1857 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1858 trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]';
1859
1860 % Build time-series test data
1861
1862 % Picks the values at random from the list
1863 fs = utils.math.randelement(fs_list, 1);
1864 nsecs = utils.math.randelement(nsecs_list, 1);
1865 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
1866 trend_0 = utils.math.randelement(trend_0_list, 1);
1867
1868 % White noise
1869 type = 'Normal';
1870 a_n1 = ao(plist('waveform', 'noise', ...
1871 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1872 a_n2 = ao(plist('waveform', 'noise', ...
1873 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1874
1875 % Constant signal
1876 a_c = ao(trend_0);
1877
1878 % Total signals
1879 a1 = a_n1 + a_c;
1880 a2 = a_n2 + a_c;
1881
1882 % Evaluate the complex coherence of the white noise time-series data
1883 win_list = specwin.getTypes;
1884 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
1885 win_type = win_type{1};
1886 switch win_type
1887 case 'Kaiser'
1888 win = specwin(win_type, 1, find(ao.getInfo('psd').plists, 'psll'));
1889 otherwise
1890 win = specwin(win_type, 1);
1891 end
1892
1893 olap = win.rov;
1894 detrend = 0;
1895 scale_type = 'C';
1896 navs = fix(utils.math.randelement(logspace(0,log10(max(a1.len/10,0)),50),1));
1897
1898 % Calculates the coherence asking for the number of averages
1899 C1 = cohere(a1, a2, plist('Win', win.type, 'olap', olap, ...
1900 'Nfft', -1, 'order', detrend, 'type', scale_type, ...
1901 'navs', navs));
1902
1903 npts_2 = find(C1.hist.plistUsed, 'Nfft');
1904 % Calculates the coherence asking for the number of points
1905 C2 = cohere(a1, a2, plist('Win', win.type, 'olap', olap, ...
1906 'Nfft', npts_2, 'order', detrend, 'type', scale_type));
1907
1908 npts_3 = fix(npts_2/2);
1909 % Calculates the coherence asking for the number of points AND the window length
1910 C3 = cohere(a1, a2, plist('Win', win.type, 'olap', olap, ...
1911 'Nfft', npts_3, ...
1912 'order', detrend, 'type', scale_type, ...
1913 'navs', navs));
1914
1915 stest = true;
1916
1917 catch err
1918 disp(err.message)
1919 stest = false;
1920 end
1921 % </SyntaxCode>
1922
1923 % <AlgoDescription>
1924 %
1925 % 1) Check that calculated objects C1 and C2 are identical
1926 % 2) Check that C3 used different values
1927 %
1928 % </AlgoDescription>
1929
1930 % <AlgoCode>
1931 atest = true;
1932
1933 if stest
1934 % Compare the navs written in the output object with the requested one
1935 if ne(C1,C2,ple3) || ...
1936 ne(find(C3.hist.plistUsed, 'Nfft'), npts_3) || eq(C3.data.navs, navs)
1937 atest = false;
1938 end
1939 else
1940 atest = false;
1941 end
1942 % </AlgoCode>
1943
1944 % Return a result structure
1945 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1946 end % END UTP_23
1947
1948 %% UTP_24
1949
1950 % <TestDescription>
1951 %
1952 % Tests that the cohere method agrees with MATLAB's mscohere when
1953 % configured to use the same parameters.
1954 %
1955 % </TestDescription>
1956 function result = utp_24
1957
1958 % <SyntaxDescription>
1959 %
1960 % Test that the applying cohere works on two AOs.
1961 %
1962 % </SyntaxDescription>
1963
1964 try
1965 % <SyntaxCode>
1966 % Construct two test AOs
1967 nsecs = 10;
1968 fs = 1000;
1969 pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');
1970 a1 = ao(pl); a2 = ao(pl);
1971 % Filter one time-series
1972 f2 = miir(plist('type', 'bandpass', 'fs', fs, 'order', 3, 'fc', [50 250]));
1973 a1f = filter(a1, plist('filter', f2));
1974 % make some cross-power
1975 a4 = a1f+a2; a4.setName;
1976 % Create the transpose of a4 to check the output data shape
1977 a4 = a4.';
1978 % Compute coherence
1979 Nfft = 2*fs;
1980 % Use different windows size as Nfft
1981 win = specwin('Hanning', 1000);
1982 pl = plist('Nfft', Nfft, 'Win', win.type, 'order', 0, 'type', 'MS');
1983 out = cohere(a4,a1,pl);
1984 % </SyntaxCode>
1985 stest = true;
1986 catch err
1987 disp(err.message)
1988 stest = false;
1989 end
1990
1991 % <AlgoDescription>
1992 %
1993 % 1) Check that output agrees with the output of MATLAB's mscohere.
1994 % 2) Check that the shape of the output data is equal to the input data
1995 %
1996 % </AlgoDescription>
1997
1998 atest = true;
1999 if stest
2000 % <AlgoCode>
2001 TOL = 1e-12;
2002
2003 % Redesign the window
2004 win = specwin('Hanning', Nfft);
2005 % Compute coherence using MATLAB's cohere
2006 [cxy, f] = mscohere(a4.y, a1.y, win.win, Nfft/2, Nfft, a1.fs);
2007 if any(abs(cxy(4:end)-out.y(4:end))>TOL), atest = false; end
2008 if ne(f, out.x), atest = false; end
2009 if ne(out, out, ple2), atest = false; end
2010 % Check the data shape
2011 if size(a4.y,1) == 1
2012 if size(out.y,1) ~= 1, atest = false; end
2013 else
2014 if size(out.y,2) ~= 1, atest = false; end
2015 end
2016 % </AlgoCode>
2017 else
2018 atest = false;
2019 end
2020
2021 % Return a result structure
2022 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2023 end % END UTP_24
2024
2025 %% UTP_25
2026
2027 % <TestDescription>
2028 %
2029 % Tests handling of units:
2030 % 1) white noise produced from normal pdf, with a given mean value and
2031 % sigma (distribution's 1st and 2nd orders)
2032 % 2) white noise produced from normal pdf, with a given mean value and
2033 % sigma (distribution's 1st and 2nd orders)
2034 % 3) complex coherence of the white noise series
2035 % 4) compares the units of the input and output
2036 %
2037
2038 % </TestDescription>
2039 function result = utp_25
2040
2041 % <SyntaxDescription>
2042 %
2043 % 1) Prepare the test tsdata:
2044 % white noise from normal distribution + offset
2045 % 2) Assign a random unit
2046 % 3) Prepare the test tsdata:
2047 % white noise from normal distribution + offset
2048 % 4) Assign a random unit
2049 % 5) complex cohere of the white noise
2050 %
2051 % </SyntaxDescription>
2052
2053 % <SyntaxCode>
2054 try
2055
2056 % Build time-series test data
2057 fs = 1;
2058 nsecs = 86400;
2059 sigma_distr_1 = 4.69e-12;
2060 mu_distr_1 = -5.11e-14;
2061 sigma_distr_2 = 6.04e-9;
2062 mu_distr_2 = 1.5e-10;
2063
2064 % White noise
2065 type = 'Normal';
2066
2067 a_n = ao(plist('waveform', 'noise', ...
2068 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1));
2069 a_const = ao(mu_distr_1);
2070 a_1 = a_n + a_const;
2071
2072 a_n = ao(plist('waveform', 'noise', ...
2073 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_2));
2074 a_const = ao(mu_distr_2);
2075 a_2 = a_n + a_const;
2076
2077 % Set units and prefix from those supported
2078 unit_list = unit.supportedUnits;
2079 % remove the first empty unit '' from the list, because then is it
2080 % possible that we add a prefix to an empty unit
2081 unit_list = unit_list(2:end);
2082 prefix_list = unit.supportedPrefixes;
2083 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
2084 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
2085
2086 % Evaluate the coherence of the time-series data
2087 win = 'Kaiser';
2088 psll = utils.math.randelement([10:10:200], 1);
2089 detrend = 0;
2090 scale_type = 'C';
2091 n_pts = nsecs*fs/10;
2092
2093 C = cohere(a_1, a_2, ...
2094 plist('Win', win, 'psll', psll, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
2095
2096 stest = true;
2097
2098 catch err
2099 disp(err.message)
2100 stest = false;
2101 end
2102 % </SyntaxCode>
2103
2104 % <AlgoDescription>
2105 %
2106 % 1) Check that (complex coherence yunits) equals [1]
2107 % 2) Check that (complex coherence xunits) equals [Hz]
2108
2109 % </AlgoDescription>
2110
2111 % <AlgoCode>
2112 atest = true;
2113
2114 if stest
2115 if ~eq(C.yunits, unit('')) || ~eq(C.xunits, unit('Hz'))
2116 atest = false;
2117 end
2118 else
2119 atest = false;
2120 end
2121 % </AlgoCode>
2122
2123 % Return a result structure
2124 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2125 end % END UTP_25
2126
2127 %% UTP_26
2128
2129 % <TestDescription>
2130 %
2131 % Tests handling of units:
2132 % 1) white noise produced from normal pdf, with a given mean value and
2133 % sigma (distribution's 1st and 2nd orders)
2134 % 2) white noise produced from normal pdf, with a given mean value and
2135 % sigma (distribution's 1st and 2nd orders)
2136 % 3) complex coherence of the white noise series
2137 % 4) compares the units of the input and output
2138 %
2139
2140 % </TestDescription>
2141 function result = utp_26
2142
2143 % <SyntaxDescription>
2144 %
2145 % 1) Prepare the test tsdata:
2146 % white noise from normal distribution + offset
2147 % 2) Assign a random unit
2148 % 3) Prepare the test tsdata:
2149 % white noise from normal distribution + offset
2150 % 4) Assign a random unit
2151 % 5) complex cohere of the white noise
2152 %
2153 % </SyntaxDescription>
2154
2155 % <SyntaxCode>
2156 try
2157
2158 % Build time-series test data
2159 fs = 1;
2160 nsecs = 86400;
2161 sigma_distr_1 = 4.69e-12;
2162 mu_distr_1 = -5.11e-14;
2163 sigma_distr_2 = 6.04e-9;
2164 mu_distr_2 = 1.5e-10;
2165
2166 % White noise
2167 type = 'Normal';
2168
2169 a_n = ao(plist('waveform', 'noise', ...
2170 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1));
2171 a_const = ao(mu_distr_1);
2172 a_1 = a_n + a_const;
2173
2174 a_n = ao(plist('waveform', 'noise', ...
2175 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_2));
2176 a_const = ao(mu_distr_2);
2177 a_2 = a_n + a_const;
2178
2179 % Set units and prefix from those supported
2180 unit_list = unit.supportedUnits;
2181 % remove the first empty unit '' from the list, because then is it
2182 % possible that we add a prefix to an empty unit
2183 unit_list = unit_list(2:end);
2184 prefix_list = unit.supportedPrefixes;
2185 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
2186 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
2187
2188 % Evaluate the coherence of the time-series data
2189 win = 'Kaiser';
2190 psll = utils.math.randelement([10:10:200], 1);
2191 detrend = 0;
2192 scale_type = 'C';
2193 n_pts = nsecs*fs/10;
2194
2195 C = cohere(a_1, a_2, ...
2196 plist('Win', win, 'psll', psll, 'Nfft', n_pts, 'order', detrend, 'type', scale_type));
2197
2198 stest = true;
2199
2200 catch err
2201 disp(err.message)
2202 stest = false;
2203 end
2204 % </SyntaxCode>
2205
2206 % <AlgoDescription>
2207 %
2208 % 1) Check that (complex coherence yunits) equals [1]
2209 % 2) Check that (complex coherence xunits) equals [Hz]
2210
2211 % </AlgoDescription>
2212
2213 % <AlgoCode>
2214 atest = true;
2215
2216 if stest
2217 if ~eq(C.yunits, unit('')) || ~eq(C.xunits, unit('Hz'))
2218 atest = false;
2219 end
2220 else
2221 atest = false;
2222 end
2223 % </AlgoCode>
2224
2225 % Return a result structure
2226 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2227 end % END UTP_26
2228
2229 %% UTP_30
2230
2231 % <TestDescription>
2232 %
2233 % Tests handling of special cases:
2234 % 1) white noise produced from normal pdf, with a given mean value and
2235 % sigma (distribution's 1st and 2nd orders)
2236 % 2) the same noise series
2237 % 3) cohere of the white noise series
2238 % 4) compares the output to unity
2239 %
2240
2241 % </TestDescription>
2242 function result = utp_30
2243
2244 % <SyntaxDescription>
2245 %
2246 % 1) Prepare the test tsdata:
2247 % white noise from normal distribution + offset
2248 % 2) Assign a random unit
2249 % 3) Prepare the test tsdata:
2250 % the same data as 1) and 2)
2251 % 4) cohere of the series
2252 %
2253 % </SyntaxDescription>
2254
2255 % <SyntaxCode>
2256 try
2257
2258 % Build time-series test data
2259 fs = 1;
2260 nsecs = 86400;
2261 sigma_distr_1 = 4.69e-12;
2262 mu_distr_1 = -5.11e-14;
2263
2264 % White noise
2265 type = 'Normal';
2266
2267 a_n = ao(plist('waveform', 'noise', ...
2268 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1));
2269 a_const = ao(mu_distr_1);
2270 a_1 = a_n + a_const;
2271
2272 % Set units and prefix from those supported
2273 unit_list = unit.supportedUnits;
2274 % remove the first empty unit '' from the list, because then is it
2275 % possible that we add a prefix to an empty unit
2276 unit_list = unit_list(2:end);
2277 prefix_list = unit.supportedPrefixes;
2278 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
2279
2280 % Build the second object as a copy of the first
2281 a_2 = a_1;
2282
2283 % Evaluate the cohere of the time-series data
2284 win = specwin('BH92');
2285 olap = win.rov;
2286 detrend = 0;
2287 n_pts = nsecs*fs/10;
2288 scale_type = 'C';
2289
2290 C = cohere(a_1, a_2, ...
2291 plist('Win', win, 'Nfft', n_pts, 'order', detrend, 'type', scale_type, 'olap', olap));
2292
2293 stest = true;
2294
2295 catch err
2296 disp(err.message)
2297 stest = false;
2298 end
2299 % </SyntaxCode>
2300
2301 % <AlgoDescription>
2302 %
2303 % 1) Check that calculated cohere equals 1
2304
2305 % </AlgoDescription>
2306
2307 % <AlgoCode>
2308 atest = true;
2309
2310 if stest
2311 if sum(ne(C.y, 1))
2312 atest = false;
2313 end
2314 else
2315 atest = false;
2316 end
2317 % </AlgoCode>
2318
2319 % Return a result structure
2320 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2321 end % END UTP_30
2322
2323 end