Mercurial > hg > ltpda

comparison testing/utp_1.1/utps/ao/utp_ao_cpsd.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_CPSD a set of UTPs for the ao/cpsd method
2 %
3 % M Hewitson 06-08-08
4 %
5 % $Id: utp_ao_cpsd.m,v 1.43 2011/07/22 11:51:46 mauro Exp $
6 %
7
8 % <MethodDescription>
9 %
10 % The cpsd method of the ao class computes the cross-spectral density between two
11 % time-series AOs.
12 %
13 % </MethodDescription>
14
15 function results = utp_ao_cpsd(varargin)
16
17 % Check the inputs
18 if nargin == 0
19
20 % Some keywords
21 class = 'ao';
22 mthd = 'cpsd';
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 against MATLAB's cpsd()
50
51 results = [results utp_11(mthd, [at1 at1], ple1)]; % Test plotinfo doesn't disappear
52
53 results = [results utp_17]; % Test units handling: CPSD
54 results = [results utp_18]; % Comparison with PSD
55 results = [results utp_24]; % Test data lengths
56 results = [results utp_25]; % Test Kaiser win and olap: CPSD
57 results = [results utp_51]; % Test number of averages: requested/obtained
58 results = [results utp_52]; % Test number of averages: correct number
59 results = [results utp_53]; % Test number of averages: syntax
60
61 disp('Done.');
62 disp('******************************************************');
63
64 elseif nargin == 1 % Check for UTP functions
65 if strcmp(varargin{1}, 'isutp')
66 results = 1;
67 else
68 results = 0;
69 end
70 else
71 error('### Incorrect inputs')
72 end
73
74 %% UTP_01
75
76 % <TestDescription>
77 %
78 % Tests that the getInfo call works for this method.
79 %
80 % </TestDescription>
81 function result = utp_01
82
83
84 % <SyntaxDescription>
85 %
86 % Test that the getInfo call works for no sets, all sets, and each set
87 % individually.
88 %
89 % </SyntaxDescription>
90
91 try
92 % <SyntaxCode>
93 % Call for no sets
94 io(1) = eval([class '.getInfo(''' mthd ''', ''None'')']);
95 % Call for all sets
96 io(2) = eval([class '.getInfo(''' mthd ''')']);
97 % Call for each set
98 for kk=1:numel(io(2).sets)
99 io(kk+2) = eval([class '.getInfo(''' mthd ''', ''' io(2).sets{kk} ''')']);
100 end
101 % </SyntaxCode>
102 stest = true;
103 catch err
104 disp(err.message)
105 stest = false;
106 end
107
108 % <AlgoDescription>
109 %
110 % 1) Check that getInfo call returned an minfo object in all cases.
111 % 2) Check that all plists have the correct parameters.
112 %
113 % </AlgoDescription>
114
115 atest = true;
116 if stest
117 % <AlgoCode>
118 % check we have minfo objects
119 if isa(io, 'minfo')
120
121 % SET 'None'
122 if ~isempty(io(1).sets), atest = false; end
123 if ~isempty(io(1).plists), atest = false; end
124 % Check all Sets
125 if ~any(strcmpi(io(2).sets, 'Default')), atest = false; end
126 if numel(io(2).plists) ~= numel(io(2).sets), atest = false; end
127 % SET 'Default'
128 if io(3).plists.nparams ~= 8, atest = false; end
129 % Check key
130 if ~io(3).plists.isparam('nfft'), atest = false; end
131 if ~io(3).plists.isparam('win'), atest = false; end
132 if ~io(3).plists.isparam('olap'), atest = false; end
133 if ~io(3).plists.isparam('order'), atest = false; end
134 if ~io(3).plists.isparam('navs'), atest = false; end
135 if ~io(3).plists.isparam('times'), atest = false; end
136 if ~io(3).plists.isparam('split'), atest = false; end
137 if ~io(3).plists.isparam('psll'), atest = false; end
138 % Check default value
139 if ~isequal(io(3).plists.find('nfft'), -1), atest = false; end
140 if ~strcmpi(io(3).plists.find('win'), defaultWinType), atest = false; end
141 if ~isequal(io(3).plists.find('olap'), -1), atest = false; end
142 if ~isequal(io(3).plists.find('order'), 0), atest = false; end
143 if ~isequal(io(3).plists.find('navs'), -1), atest = false; end
144 if ~isEmptyDouble(io(3).plists.find('times')), atest = false; end
145 if ~isEmptyDouble(io(3).plists.find('split')), atest = false; end
146 if ~isequal(io(3).plists.find('psll'), 200), atest = false; end
147 % Check options
148 if ~isequal(io(3).plists.getOptionsForParam('nfft'), {-1}), atest = false; disp('1'); end
149 if ~isequal(io(3).plists.getOptionsForParam('win'), specwin.getTypes), atest = false;disp('2'); end
150 if ~isequal(io(3).plists.getOptionsForParam('olap'), {-1}), atest = false; disp('3');end
151 if ~isequal(io(3).plists.getOptionsForParam('order'), {-1 0 1 2 3 4 5 6 7 8 9}), atest = false;disp('4'); end
152 if ~isequal(io(3).plists.getOptionsForParam('navs'), {-1}), atest = false;disp('5'); end
153 if ~isequal(io(3).plists.getOptionsForParam('times'), {[]}), atest = false;disp('6'); end
154 if ~isequal(io(3).plists.getOptionsForParam('split'), {[]}), atest = false;disp('6'); end
155 if ~isequal(io(3).plists.getOptionsForParam('psll'), {200}), atest = false;disp('7'); end
156 end
157 % </AlgoCode>
158 else
159 atest = false;
160 end
161
162 % Return a result structure
163 result = utp_prepare_result(atest, stest, dbstack, mfilename);
164 end % END UTP_01
165
166 %% UTP_02
167
168 % <TestDescription>
169 %
170 % Tests that the cpsd method works with a vector of AOs as input. (only
171 % with two objects in the vector)
172 %
173 % </TestDescription>
174 function result = utp_02
175
176 % <SyntaxDescription>
177 %
178 % Test that the cpsd method works for a vector of AOs as input.
179 %
180 % </SyntaxDescription>
181
182 try
183 % <SyntaxCode>
184 avec = [at1 at5];
185 out = cpsd(avec);
186 % </SyntaxCode>
187 stest = true;
188 catch err
189 disp(err.message)
190 stest = false;
191 end
192
193 % <AlgoDescription>
194 %
195 % 1) Check that the number of elements in 'out' is equal to 1
196 %
197 % </AlgoDescription>
198
199 atest = true;
200 if stest
201 % <AlgoCode>
202 % Check we have the correct number of outputs
203 if numel(out) ~= 1, atest = false; end
204 % </AlgoCode>
205 else
206 atest = false;
207 end
208
209 % Return a result structure
210 result = utp_prepare_result(atest, stest, dbstack, mfilename);
211 end % END UTP_02
212
213 %% UTP_03
214
215 % <TestDescription>
216 %
217 % Tests that the cpsd method doesn't work with a matrix of AOs as input.
218 %
219 % </TestDescription>
220 function result = utp_03
221
222 % <SyntaxDescription>
223 %
224 % Test that the cpsd method doesn't work for a matrix of AOs as input.
225 %
226 % </SyntaxDescription>
227
228 try
229 % <SyntaxCode>
230 amat = [at1 at5 at6; at5 at6 at1];
231 out = cpsd(amat);
232 % </SyntaxCode>
233 stest = false;
234 catch err
235 stest = true;
236 end
237
238 % <AlgoDescription>
239 %
240 % 1) Nothing to check.
241 %
242 % </AlgoDescription>
243
244 atest = true;
245 if stest
246 % <AlgoCode>
247 % </AlgoCode>
248 else
249 atest = false;
250 end
251
252 % Return a result structure
253 result = utp_prepare_result(atest, stest, dbstack, mfilename);
254 end % END UTP_03
255
256 %% UTP_04
257
258 % <TestDescription>
259 %
260 % Tests that the cpsd method works with a list of AOs as input.
261 %
262 % </TestDescription>
263 function result = utp_04
264
265 % <SyntaxDescription>
266 %
267 % Test that the cpsd method works for a list of AOs as input.
268 %
269 % </SyntaxDescription>
270
271 try
272 % <SyntaxCode>
273 out = cpsd(at1,at5);
274 % </SyntaxCode>
275 stest = true;
276 catch err
277 disp(err.message)
278 stest = false;
279 end
280
281 % <AlgoDescription>
282 %
283 % 1) Check that the number of elements in 'out' is equal to 1
284 %
285 % </AlgoDescription>
286
287 atest = true;
288 if stest
289 % <AlgoCode>
290 % Check we have the correct number of outputs
291 if numel(out) ~= 1, atest = false; end
292 % </AlgoCode>
293 else
294 atest = false;
295 end
296
297 % Return a result structure
298 result = utp_prepare_result(atest, stest, dbstack, mfilename);
299 end % END UTP_04
300
301 %% UTP_05
302
303 % <TestDescription>
304 %
305 % Tests that the cpsd method doesn't work with a mix of different shaped
306 % AOs as input.
307 %
308 % </TestDescription>
309 function result = utp_05
310
311 % <SyntaxDescription>
312 %
313 % Test that the cpsd method doesn't work with an input of matrices and
314 % vectors and single AOs.
315 %
316 % </SyntaxDescription>
317
318 try
319 % <SyntaxCode>
320 out = cpsd([at5 at6],[at5 at1; at6 at1],at6);
321 stest = false;
322 % </SyntaxCode>
323 catch err
324 stest = true;
325 end
326
327 % <AlgoDescription>
328 %
329 % 1) Nothing to check
330 %
331 % </AlgoDescription>
332
333 atest = true;
334 if stest
335 % <AlgoCode>
336 % </AlgoCode>
337 else
338 atest = false;
339 end
340
341 % Return a result structure
342 result = utp_prepare_result(atest, stest, dbstack, mfilename);
343 end % END UTP_05
344
345 %% UTP_06
346
347 % <TestDescription>
348 %
349 % Tests that the cpsd method properly applies history.
350 %
351 % </TestDescription>
352 function result = utp_06
353
354 % <SyntaxDescription>
355 %
356 % Test that the result of applying the cpsd method can be processed back
357 % to an m-file.
358 %
359 % </SyntaxDescription>
360
361 try
362 % <SyntaxCode>
363 out = cpsd(at5,at6);
364 mout = rebuild(out);
365 % </SyntaxCode>
366 stest = true;
367 catch err
368 disp(err.message)
369 stest = false;
370 end
371
372 % <AlgoDescription>
373 %
374 % 1) Check that the last entry in the history of 'out' corresponds to
375 % 'cpsd'.
376 % 2) Check that the re-built object is the same as 'out'.
377 %
378 % </AlgoDescription>
379
380 atest = true;
381 if stest
382 % <AlgoCode>
383 % Check the last step in the history of 'out'
384 if ~strcmp(out.hist.methodInfo.mname, 'cpsd'), atest = false; end
385 % Check the re-built object
386 if ~eq(mout, out, ple2), atest = false; end
387 % </AlgoCode>
388 else
389 atest = false;
390 end
391
392 % Return a result structure
393 result = utp_prepare_result(atest, stest, dbstack, mfilename);
394 end % END UTP_06
395
396 %% UTP_07
397
398 % <TestDescription>
399 %
400 % Tests that the cpsd method can not modify the input AO.
401 %
402 % </TestDescription>
403 function result = utp_07
404
405 % <SyntaxDescription>
406 %
407 % Test that the cpsd method can not modify the input AO.
408 % The method must throw an error for the modifier call.
409 %
410 % </SyntaxDescription>
411
412 try
413 % <SyntaxCode>
414 % copy at1 to work with
415 ain = ao(at1);
416 % modify ain
417 ain.cpsd(at5);
418 % </SyntaxCode>
419 stest = false;
420 catch err
421 stest = true;
422 end
423
424 % <AlgoDescription>
425 %
426 % 1) Nothing to check.
427 %
428 % </AlgoDescription>
429
430 atest = true;
431 if stest
432 % <AlgoCode>
433 % </AlgoCode>
434 else
435 atest = false;
436 end
437
438 % Return a result structure
439 result = utp_prepare_result(atest, stest, dbstack, mfilename);
440 end % END UTP_07
441
442 %% UTP_08
443
444 % <TestDescription>
445 %
446 % Test the shape of the output.
447 %
448 % </TestDescription>
449 function result = utp_08
450
451 % <SyntaxDescription>
452 %
453 % Test that the cpsd method keeps the data shape of the input object. The
454 % input AO must be an AO with row data and an AO with column data.
455 %
456 % </SyntaxDescription>
457
458 try
459 % <SyntaxCode>
460 out1 = cpsd(at5, at6);
461 out2 = cpsd(at6, at5);
462 % </SyntaxCode>
463 stest = true;
464 catch err
465 disp(err.message)
466 stest = false;
467 end
468
469 % <AlgoDescription>
470 %
471 % 1) Check that the shape of the output data doesn't change.
472 %
473 % </AlgoDescription>
474
475 atest = true;
476 if stest
477 % <AlgoCode>
478 % Check the shape of the output data
479 if size(out1.data.y, 2) ~= 1, atest = false; end
480 if size(out2.data.y, 1) ~= 1, atest = false; end
481 % </AlgoCode>
482 else
483 atest = false;
484 end
485
486 % Return a result structure
487 result = utp_prepare_result(atest, stest, dbstack, mfilename);
488 end % END UTP_08
489
490 %% UTP_09
491
492 % <TestDescription>
493 %
494 % Check that the cpsd method pass back the output objects to a list of
495 % output variables or to a single variable.
496 %
497 % </TestDescription>
498 function result = utp_09
499
500 % <SyntaxDescription>
501 %
502 % This test is not longer necessary because the cpsd method pass back
503 % always only one object.
504 %
505 % </SyntaxDescription>
506
507 try
508 % <SyntaxCode>
509 % </SyntaxCode>
510 stest = true;
511 catch err
512 disp(err.message)
513 stest = false;
514 end
515
516 % <AlgoDescription>
517 %
518 % 1) Nothing to check.
519 %
520 % </AlgoDescription>
521
522 atest = true;
523 if stest
524 % <AlgoCode>
525 % </AlgoCode>
526 else
527 atest = false;
528 end
529
530 % Return a result structure
531 result = utp_prepare_result(atest, stest, dbstack, mfilename);
532 end % END UTP_09
533
534 %% UTP_10
535
536 % <TestDescription>
537 %
538 % Tests that the cpsd method agrees with MATLAB's cpsd when
539 % configured to use the same parameters.
540 %
541 % </TestDescription>
542 function result = utp_10
543
544 % <SyntaxDescription>
545 %
546 % Test that applying cpsd works on two AOs.
547 %
548 % </SyntaxDescription>
549
550 try
551 % <SyntaxCode>
552 % Construct two test AOs
553 nsecs = 10;
554 fs = 1000;
555 pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');
556 a1 = ao(pl); a2 = ao(pl);
557 % Filter one time-series
558 f2 = miir(plist('type', 'bandpass', 'fs', fs, 'order', 3, 'fc', [50 250]));
559 a1f = filter(a1, plist('filter', f2));
560 % make some cross-power
561 a4 = a1f+a2; a4.setName;
562 % Compute cpsd
563 Nfft = 2*fs;
564 win = specwin('Hanning', Nfft);
565 pl = plist('Nfft', Nfft, 'Win', win.type, 'order', -1);
566 out = cpsd(a4,a1,pl);
567 % </SyntaxCode>
568 stest = true;
569 catch err
570 disp(err.message)
571 stest = false;
572 end
573
574 % <AlgoDescription>
575 %
576 % 1) Check that output agrees with the output of MATLAB's cpsd.
577 %
578 % </AlgoDescription>
579
580 atest = true;
581 if stest
582 % <AlgoCode>
583 % Compute cpsd using MATLAB's cpsd
584 [cxy, f] = cpsd(a4.data.y, a1.data.y, win.win, Nfft/2, Nfft, a1.data.fs);
585 if ~utils.math.isequal(cxy(:), out.data.y(:)) || ~utils.math.isequal(f, out.data.getX), atest = false; end
586 % </AlgoCode>
587 else
588 atest = false;
589 end
590
591 % Return a result structure
592 result = utp_prepare_result(atest, stest, dbstack, mfilename);
593 end % END UTP_10
594
595
596 %% UTP_17
597
598 % <TestDescription>
599 %
600 % Tests handling of units:
601 % 1) white noise produced from normal pdf, with a given mean value and
602 % sigma (distribution's 1st and 2nd orders)
603 % 2) white noise produced from normal pdf, with a given mean value and
604 % sigma (distribution's 1st and 2nd orders)
605 % 3) CPSD of the white noise series
606 % 4) compares the units of the input and output
607 %
608
609 % </TestDescription>
610 function result = utp_17
611
612 % <SyntaxDescription>
613 %
614 % 1) Prepare the test tsdata:
615 % white noise from normal distribution + offset
616 % 2) Assign a random unit
617 % 3) Prepare the test tsdata:
618 % white noise from normal distribution + offset
619 % 4) Assign a random unit
620 % 5) CPSD of the white noise
621 %
622 % </SyntaxDescription>
623
624 % <SyntaxCode>
625 try
626
627 noise_type = 'Normal';
628 win_type = 'BH92';
629
630 [a_1, a_2, spec, spec1] = prepare_analyze_noise(win_type, noise_type, plist);
631
632 stest = true;
633
634 catch err
635 disp(err.message)
636 stest = false;
637 end
638 % </SyntaxCode>
639
640 % <AlgoDescription>
641 %
642 % 1) Check that (calculated CPSD yunits) equals
643 % input_1 units*input_2 units/Hz
644 %
645 % </AlgoDescription>
646
647 % <AlgoCode>
648 atest = true;
649 u = simplifyYunits(a_1.* a_2, plist('prefixes', false, 'exceptions', 'Hz'));
650 if stest
651 if ne(spec.Cxy.yunits, u.yunits * unit('Hz^-1')) || ne(spec.Cxy.xunits, unit('Hz'))
652 atest = false;
653 end
654 if ne(spec.Cyx.yunits, u.yunits * unit('Hz^-1')) || ne(spec.Cyx.xunits, unit('Hz'))
655 atest = false;
656 end
657 else
658 atest = false;
659 end
660 % </AlgoCode>
661
662 % Return a result structure
663 result = utp_prepare_result(atest, stest, dbstack, mfilename);
664 end % END UTP_17
665
666 %% UTP_18
667
668 % <TestDescription>
669 %
670 % Tests handling of units:
671 % 1) white noise produced from normal pdf, with a given mean value and
672 % sigma (distribution's 1st and 2nd orders)
673 % 2) white noise produced from normal pdf, with a given mean value and
674 % sigma (distribution's 1st and 2nd orders)
675 % 3) CPSD of the white noise series
676 %
677 % Comparison with PSD:
678 % 4) compares the off-diagonal terms to check they are complex-conjugated
679 % 5) compares the diagonal terms with PSD of the individual noise
680 %
681
682 % </TestDescription>
683 function result = utp_18
684
685 % <SyntaxDescription>
686 %
687 % 1) Prepare the test tsdata:
688 % white noise from normal distribution + offset
689 % 2) Assign a random unit
690 % 3) Prepare the test tsdata:
691 % white noise from normal distribution + offset
692 % 4) Assign a random unit
693 % 5) CPSD of the white noise
694 % 6) PSD of the white noise
695 %
696 % </SyntaxDescription>
697
698 % <SyntaxCode>
699 try
700
701 noise_type = 'Uniform';
702 win_type = 'BH92';
703
704 [a_1, a_2, spec, spec2] = prepare_analyze_noise(win_type, noise_type, plist);
705
706 stest = true;
707
708 catch err
709 disp(err.message)
710 stest = false;
711 end
712 % </SyntaxCode>
713
714 % <AlgoDescription>
715 %
716 % 1) Check that CPSD(x,y) equals conj(CPSD(y,x))
717 % 2) Check that CPSD(x,x) equals PSD(x)
718 % 3) Check that CPSD(y,y) equals PSD(y)
719 %
720 % </AlgoDescription>
721
722 % <AlgoCode>
723 atest = true;
724
725 if stest
726 if ne(spec.Cxy.y, conj(spec.Cyx.y)), atest = false; end
727 if ne(spec.Cxy.x, spec.Cyx.x), atest = false; end
728 if ne(spec.Cxx.data, spec.S_1.data), atest = false; end
729 if ne(spec.Cyy.data, spec.S_2.data), atest = false; end
730 else
731 atest = false;
732 end
733 % </AlgoCode>
734
735 % Return a result structure
736 result = utp_prepare_result(atest, stest, dbstack, mfilename);
737 end % END UTP_18
738
739
740 %% UTP_24
741
742 % <TestDescription>
743 %
744 % Tests that differently sized data sets are treated properly
745 %
746 % </TestDescription>
747 function result = utp_24
748
749 % <SyntaxDescription>
750 %
751 % Test that applying cpsd works on two AOs.
752 %
753 % </SyntaxDescription>
754
755 try
756 % <SyntaxCode>
757 % Construct two test AOs
758 nsecs = [10000:1:20000];
759 fs = 1;
760 pl = plist('fs', fs, 'tsfcn', 'randn(size(t))');
761 a1 = ao(pl.pset('nsecs', utils.math.randelement(nsecs, 1)));
762 a2 = ao(pl.pset('nsecs', utils.math.randelement(nsecs, 1)));
763 len_1 = a1.len;
764 len_2 = a2.len;
765 % Filter one time-series
766 f2 = miir(plist('type', 'bandpass', 'fs', fs, 'order', 3, 'fc', [.050 .25]));
767 a1f = filter(a1, plist('filter', f2));
768 % Compute cpsd
769 Nfft = -1;
770 win = 'Hanning';
771 pl = plist('Nfft', Nfft, 'Win', win, 'order', -1);
772 out = cpsd(a2,a1f,pl);
773 % </SyntaxCode>
774 stest = true;
775 catch err
776 disp(err.message)
777 stest = false;
778 end
779
780 % <AlgoDescription>
781 %
782 % 1) Check that cpsd used the length of the shortest ao.
783 %
784 % </AlgoDescription>
785
786 atest = true;
787 if stest
788 % <AlgoCode>
789 % Compare the nfft with the length of the input data
790
791 if out.x(2) ~= 1/min(len_1,len_2)
792 atest = false;
793 end
794 % </AlgoCode>
795 else
796 atest = false;
797 end
798
799 % Return a result structure
800 result = utp_prepare_result(atest, stest, dbstack, mfilename);
801 end % END UTP_24
802
803 %% UTP_25
804
805 % <TestDescription>
806 %
807 % Tests handling of units:
808 % 1) white noise produced from normal pdf, with a given mean value and
809 % sigma (distribution's 1st and 2nd orders)
810 % 2) white noise produced from normal pdf, with a given mean value and
811 % sigma (distribution's 1st and 2nd orders)
812 % 3) CPSD of the white noise series
813 % 4) compares the units of the input and output
814 %
815
816 % </TestDescription>
817 function result = utp_25
818
819 % <SyntaxDescription>
820 %
821 % 1) Prepare the test tsdata:
822 % white noise from normal distribution + offset
823 % 2) Assign a random unit
824 % 3) Prepare the test tsdata:
825 % white noise from normal distribution + offset
826 % 4) Assign a random unit
827 % 5) CPSD of the white noise
828 %
829 % </SyntaxDescription>
830
831 % <SyntaxCode>
832 try
833
834 % Build time-series test data
835 fs = 1;
836 nsecs = 86400;
837 sigma_distr_1 = 4.69e-12;
838 mu_distr_1 = -5.11e-14;
839 sigma_distr_2 = 6.04e-9;
840 mu_distr_2 = 1.5e-10;
841
842 % White noise
843 type = 'Normal';
844
845 a_n = ao(plist('waveform', 'noise', ...
846 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1));
847 a_const = ao(mu_distr_1);
848 a_1 = a_n + a_const;
849
850 a_n = ao(plist('waveform', 'noise', ...
851 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_2));
852 a_const = ao(mu_distr_2);
853 a_2 = a_n + a_const;
854
855 % Set units and prefix from those supported
856 unit_list = unit.supportedUnits;
857 % remove the first empty unit '' from the list, because then is it
858 % possible that we add a prefix to an empty unit
859 unit_list = unit_list(2:end);
860 prefix_list = unit.supportedPrefixes;
861 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
862 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
863
864 % Evaluate the cpsd of the time-series data, using Kaiser window
865 % Psll and olap are not set
866 win = ('Kaiser');
867 detrend = 0;
868 n_pts = nsecs*fs/10;
869
870 C = cpsd(a_1, a_2, plist('Win', win, 'Nfft', n_pts, 'order', detrend));
871
872 stest = true;
873
874 catch err
875 disp(err.message)
876 stest = false;
877 end
878 % </SyntaxCode>
879
880 % <AlgoDescription>
881 %
882 % 1) Check that (calculated CPSD yunits) equals
883 %input_1 units*input_2 units/Hz
884
885 % </AlgoDescription>
886
887 % <AlgoCode>
888 atest = true;
889 u = simplifyYunits(a_1.* a_2, plist('prefixes', false, 'exceptions', 'Hz'));
890 if stest
891 if ne(C.yunits, u.yunits * unit('Hz^-1')) || ne(C.xunits, unit('Hz'))
892 atest = false;
893 end
894 else
895 atest = false;
896 end
897 % </AlgoCode>
898
899 % Return a result structure
900 result = utp_prepare_result(atest, stest, dbstack, mfilename);
901 end % END UTP_25
902
903 %% UTP_51
904
905 % <TestDescription>
906 %
907 % Tests the possibility to set the number of averages rather than setting the Nfft:
908 % 1) white noise produced from normal pdf, with:
909 % a given mean value and sigma (distribution's 1st and 2nd order)
910 % 2) cpsd of the noise, without detrending, random window, set number of
911 % averages
912 % 3) check the effective number of averages
913 %
914
915 % </TestDescription>
916 function result = utp_51
917
918 % <SyntaxDescription>
919 %
920 % 1) Prepare the test tsdata:
921 % white noise from normal distribution + offset
922 % 2) cpsd of the noise, without detrending, random window, set number of
923 % averages
924 %
925 % </SyntaxDescription>
926
927 % <SyntaxCode>
928 try
929
930 noise_type = 'Normal';
931
932 % Evaluate the cpsd of the white noise time-series data
933 [a_1, a_2, C1, C2, navs] = prepare_analyze_noise_navs(noise_type, plist);
934
935 stest = true;
936
937 catch err
938 disp(err.message)
939 stest = false;
940 end
941 % </SyntaxCode>
942
943 % <AlgoDescription>
944 %
945 % 1) Check that calculated navs are identical to those requested
946 %
947 % </AlgoDescription>
948
949 % <AlgoCode>
950 atest = true;
951
952 if stest
953 if ne(navs, C1.data.navs)
954 if ne(find(C1.hist.plistUsed, 'navs'), C1.data.navs)
955 atest = false;
956 end
957 end
958 else
959 atest = false;
960 end
961 % </AlgoCode>
962
963 % Return a result structure
964 result = utp_prepare_result(atest, stest, dbstack, mfilename);
965 end % END UTP_51
966
967 %% UTP_52
968
969 % <TestDescription>
970 %
971 % Tests the possibility to set the number of averages rather than setting the Nfft:
972 % 1) white noise produced from uniform pdf, with:
973 % a given mean value and sigma (distribution's 1st and 2nd order)
974 % 2) cpsd of the noise, without detrending, random window, random navs
975 % 3) get the number of averages
976 % 4) get the nfft used
977 % 5) run cpsd again, with the nfft used
978 % 6) compare the calculated objects
979 %
980
981 % </TestDescription>
982 function result = utp_52
983
984 % <SyntaxDescription>
985 %
986 % 1) white noise produced from uniform pdf, with:
987 % a given mean value and sigma (distribution's 1st and 2nd order)
988 % 2) cpsd of the noise, without detrending, random window, random navs
989 % 3) get the number of averages
990 % 4) get the nfft used
991 % 5) run cpsd again, with the nfft used
992 %
993 % </SyntaxDescription>
994
995 % <SyntaxCode>
996 try
997
998 noise_type = 'Uniform';
999
1000 % Evaluate the cpsd of the white noise time-series data
1001 [a_1, a_2, C1, C2, navs] = prepare_analyze_noise_navs(noise_type, plist);
1002
1003 stest = true;
1004
1005 catch err
1006 disp(err.message)
1007 stest = false;
1008 end
1009 % </SyntaxCode>
1010
1011 % <AlgoDescription>
1012 %
1013 % 1) Check that calculated objects C1 and C2 are identical
1014 %
1015 % </AlgoDescription>
1016
1017 % <AlgoCode>
1018 atest = true;
1019
1020 if stest
1021 % Compare the output objects
1022 if ne(C1, C2, ple3)
1023 atest = false;
1024 end
1025 else
1026 atest = false;
1027 end
1028 % </AlgoCode>
1029
1030 % Return a result structure
1031 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1032 end % END UTP_52
1033
1034 %% UTP_53
1035
1036 % <TestDescription>
1037 %
1038 % Tests the possibility to set the number of averages rather than setting the Nfft:
1039 % 1) white noise produced from normal pdf, with:
1040 % a given mean value and sigma (distribution's 1st and 2nd order)
1041 % 2) cpsd of the noise, without detrending, random window, random navs
1042 % 3) get the number of averages
1043 % 4) get the nfft used
1044 % 5) run cpsd again, with the nfft used
1045 % 6) compare navs, nfft, psds
1046 %
1047
1048 % </TestDescription>
1049 function result = utp_53
1050
1051 % <SyntaxDescription>
1052 %
1053 % 1) white noise produced from normal pdf, with:
1054 % a given mean value and sigma (distribution's 1st and 2nd order)
1055 % 2) cpsd of the noise, without detrending, random window, random navs
1056 % 3) get the number of averages
1057 % 4) get the nfft used
1058 % 5) run cpsd again, with the nfft used
1059 % 6) run cpsd again, with conflicting parameters, and verify it uses
1060 % nfft rather than navs
1061 %
1062 % </SyntaxDescription>
1063
1064 % <SyntaxCode>
1065 try
1066
1067 noise_type = 'Uniform';
1068
1069 % Evaluate the cpsd of the white noise time-series data
1070 [a_1, a_2, C1, C2, navs] = prepare_analyze_noise_navs(noise_type, plist);
1071
1072 npts_3 = fix(find(C1.hist.plistUsed, 'Nfft')/2);
1073
1074 % Calculates the cpsd asking for the number of points AND the window length
1075 pl_spec = C1.hist.plistUsed;
1076 pl_spec.pset('Nfft', npts_3, 'navs', navs);
1077 C3 = cpsd(a_1, a_2, pl_spec);
1078
1079 stest = true;
1080
1081 catch err
1082 disp(err.message)
1083 stest = false;
1084 end
1085 % </SyntaxCode>
1086
1087 % <AlgoDescription>
1088 %
1089 % 1) Check that calculated objects C1 and C2 are identical
1090 % 2) Check that C3 used different values
1091 %
1092 % </AlgoDescription>
1093
1094 % <AlgoCode>
1095 atest = true;
1096
1097 if stest
1098 % Compare the navs written in the output object with the requested one
1099 if ne(C1,C2,ple3) || ...
1100 ne(find(C3.hist.plistUsed, 'Nfft'), npts_3) || eq(C3.data.navs, navs)
1101 atest = false;
1102 end
1103 else
1104 atest = false;
1105 end
1106 % </AlgoCode>
1107
1108 % Return a result structure
1109 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1110 end % END UTP_53
1111
1112 %% Helper function for window call construction
1113
1114 function [a_1, a_2, spec1, spec2] = prepare_analyze_noise(win_type, noise_type, pli)
1115 % Array of parameters to pick from
1116 fs_list = [0.1;1;2;5;10];
1117 nsecs_list = [20 100 1000:1000:10000]';
1118 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1119 trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]';
1120
1121
1122 % Build time-series test data
1123
1124 % Picks the values at random from the list
1125 fs = utils.math.randelement(fs_list, 1);
1126 nsecs = utils.math.randelement(nsecs_list, 1);
1127 sigma_distr_1 = utils.math.randelement(sigma_distr_list, 1);
1128 sigma_distr_2 = utils.math.randelement(sigma_distr_list, 1);
1129 trend_0_1 = utils.math.randelement(trend_0_list, 1);
1130 trend_0_2 = utils.math.randelement(trend_0_list, 1);
1131
1132 % Pick units and prefix from those supported
1133 unit_list = unit.supportedUnits;
1134 % remove the first empty unit '' from the list, because then is it
1135 % possible that we add a prefix to an empty unit
1136 unit_list = unit_list(2:end);
1137 prefix_list = unit.supportedPrefixes;
1138
1139 % White noise
1140 a_n = ao(plist('waveform', 'noise', ...
1141 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_1));
1142
1143 % Constant signal
1144 a_c = ao(trend_0_1);
1145
1146 % Total signal
1147 a_1 = a_n + a_c;
1148
1149 % White noise
1150 a_n = ao(plist('waveform', 'noise', ...
1151 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr_2));
1152 % Constant signal
1153 a_c = ao(trend_0_2);
1154
1155 % Total signal
1156 a_2 = a_n + a_c;
1157
1158 % Set units
1159 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1160 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1161
1162 % Evaluate the cpsd of the white noise time-series data
1163 olap = 0;
1164 detrend_order = 0;
1165
1166 switch lower(win_type)
1167 case 'kaiser'
1168 psll = find(pli, 'psll');
1169 if isempty(psll)
1170 psll = find(ao.getInfo('psd').plists, 'psll');
1171 end
1172 pl_spec = plist('Win', win_type, 'psll', psll, 'olap', olap, 'order', detrend_order);
1173
1174 case 'levelledhanning'
1175 levelCoef = find(pli, 'levelCoef');
1176 if isempty(levelCoef)
1177 levelCoef = 1;
1178 end
1179 pl_spec = plist('Win', win_type, 'levelCoef', levelCoef, 'olap', olap, 'order', detrend_order);
1180
1181 otherwise
1182 pl_spec = plist('Win', win_type, 'olap', olap, 'order', detrend_order);
1183
1184 end
1185
1186 if find(pli, 'win_obj')
1187 % Calls the cpsd applying the detrend and window internally
1188 % (passig window object)
1189 spec2.pl = pl_spec;
1190 spec2.Cxy = cpsd(a_1, a_2, spec2.pl);
1191 spec2.Cyx = cpsd(a_2, a_1, spec2.pl);
1192 spec2.Cxx = cpsd(a_1, a_1, spec2.pl);
1193 spec2.Cyy = cpsd(a_2, a_2, spec2.pl);
1194 spec2.S_1 = simplifyYunits(psd(a_1, spec2.pl), ...
1195 plist('prefixes', false, 'exceptions','Hz'));
1196 spec2.S_2 = simplifyYunits(psd(a_2, spec2.pl), ...
1197 plist('prefixes', false, 'exceptions','Hz'));
1198 else
1199 spec2 = struct;
1200 end
1201 % Calls the cpsd applying the detrend and window internally
1202 % (passig window name)
1203 spec1.pl = pl_spec.pset('Win', win_type);
1204 spec1.Cxy = cpsd(a_1, a_2, spec1.pl);
1205 spec1.Cyx = cpsd(a_2, a_1, spec1.pl);
1206 spec1.Cxx = cpsd(a_1, a_1, spec1.pl);
1207 spec1.Cyy = cpsd(a_2, a_2, spec1.pl);
1208 spec1.S_1 = simplifyYunits(psd(a_1, spec1.pl), ...
1209 plist('prefixes', false, 'exceptions','Hz'));
1210 spec1.S_2 = simplifyYunits(psd(a_2, spec1.pl), ...
1211 plist('prefixes', false, 'exceptions','Hz'));
1212
1213 end
1214
1215 %% Helper function for window call construction, navs option
1216
1217 function [a_1, a_2, C1, C2, navs] = prepare_analyze_noise_navs(noise_type, pli)
1218 % Array of parameters to pick from
1219 fs_list = [0.1;1;2;5;10];
1220 nsecs_list = [2000:1000:10000]';
1221 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
1222 trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]';
1223
1224 % Build time-series test data
1225
1226 % Picks the values at random from the list
1227 fs = utils.math.randelement(fs_list, 1);
1228 nsecs = utils.math.randelement(nsecs_list, 1);
1229 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
1230 trend_0 = utils.math.randelement(trend_0_list, 1);
1231
1232 % Pick units and prefix from those supported
1233 unit_list = unit.supportedUnits;
1234 % remove the first empty unit '' from the list, because then is it
1235 % possible that we add a prefix to an empty unit
1236 unit_list = unit_list(2:end);
1237 prefix_list = unit.supportedPrefixes;
1238
1239 % White noise
1240 a_n1 = ao(plist('waveform', 'noise', ...
1241 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1242 a_n2 = ao(plist('waveform', 'noise', ...
1243 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
1244
1245 % Constant signal
1246 a_c = ao(trend_0);
1247
1248 % Total signals
1249 a_1 = a_n1 + a_c;
1250 a_2 = a_n2 + a_c;
1251
1252 % Set units
1253 a_1.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1254 a_2.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
1255
1256 % Evaluate the cpsd of the white noise time-series data
1257 olap = 0;
1258 detrend_order = 0;
1259 n_pts = -1;
1260
1261 navs = fix(utils.math.randelement(logspace(0,log10(max(0,a_1.len/10)),50),1));
1262
1263 % Evaluate the cpsd of the white noise time-series data
1264 % Window
1265 win_list = specwin.getTypes;
1266 win_type = utils.math.randelement(win_list(~strcmpi(win_list, 'levelledhanning')), 1);
1267 win_type = win_type{1};
1268
1269 switch lower(win_type)
1270 case 'kaiser'
1271 psll = utils.math.randelement([0:10:200],1);
1272 if psll == 0
1273 psll = find(ao.getInfo('psd').plists, 'psll');
1274 end
1275 pl_spec = plist('Win', win_type, 'psll', psll, 'olap', olap, 'order', detrend_order);
1276 otherwise
1277 pl_spec = plist('Win', win_type, 'olap', olap, 'order', detrend_order);
1278 end
1279
1280 % Calls cpsd asking for the number of averages
1281 pl_spec.pset('Nfft', n_pts, 'navs', navs);
1282 C1 = cpsd(a_1, a_2, pl_spec);
1283
1284 % Calls cpsd asking for the number of points just evaluated
1285 pl_spec.pset('Nfft', find(C1.hist.plistUsed, 'Nfft'));
1286 pl_spec.remove('navs');
1287 C2 = cpsd(a_1, a_2, pl_spec);
1288
1289 end
1290
1291 end