Mercurial > hg > ltpda

comparison testing/utp_1.1/utps/ao/utp_ao_psd.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_PSD a set of UTPs for the ao/psd method
2 %
3 % M Hewitson 06-08-08
4 %
5 % $Id: utp_ao_psd.m,v 1.48 2011/07/21 06:46:48 mauro Exp $
6 %
7
8 % <MethodDescription>
9 %
10 % The psd method of the ao class computes the spectral density of time-series AOs.
11 %
12 % </MethodDescription>
13
14 function results = utp_ao_psd(varargin)
15
16 % Check the inputs
17 if nargin == 0
18
19 addpath(fullfile(fileparts(which(mfilename)), 'reference_files'))
20
21 % Some keywords
22 class = 'ao';
23 mthd = 'psd';
24
25 results = [];
26 disp('******************************************************');
27 disp(['**** Running UTPs for ' class '/' mthd]);
28 disp('******************************************************');
29
30 % Test AOs
31 [at1,at2,at3,at4,at5,at6,atvec,atmat] = eval(['get_test_objects_' class]);
32
33 % Exception list for the UTPs:
34 [ple1,ple2,ple3,ple4,ple5,ple6] = get_test_ples();
35
36 % Get default window from the preferences
37 prefs = getappdata(0, 'LTPDApreferences');
38 defaultWinType = char(prefs.getMiscPrefs.getDefaultWindow);
39
40 % Very-long white noise reference time-series
41 a_n_long = ao(plist('filename', 'ref_whitenoise_20000s_10Hz.xml'));
42
43 % Split the reference data into different segments
44 nsecs = a_n_long.data.nsecs;
45 Nsegments = 100;
46
47 ap = a_n_long.split(plist('split_type', 'chunks', 'chunks', Nsegments));
48
49 % Run the tests
50 results = [results utp_01]; % getInfo call
51 results = [results utp_02]; % Vector input
52 results = [results utp_03]; % Matrix input
53 results = [results utp_04]; % List input
54 results = [results utp_05]; % Test with mixed input
55 results = [results utp_06]; % Test history is working
56 results = [results utp_07]; % Test against MATLAB's pwelch
57 results = [results utp_08]; % Test output of the data
58 results = [results utp_09]; % Test against MATLAB's pwelch
59 results = [results utp_10]; % Test against MATLAB's pwelch
60
61 results = [results utp_11(mthd, at1, ple1)]; % Test plotinfo doesn't disappear
62
63 results = [results utp_12]; % Test "conservation of energy":
64 % - white noise from uniform pdf, random parameters
65 % - no detrending, Rectangular window
66 results = [results utp_13]; % Test "conservation of energy":
67 % - white noise from normal pdf, fixed parameters
68 % - no detrending, Rectangular window
69 results = [results utp_14]; % Test "conservation of energy":
70 % - white noise from uniform pdf, fixed parameters
71 % - no detrending, Rectangular window
72 results = [results utp_15]; % Test "conservation of energy":
73 % - white noise from normal pdf, fixed parameters
74 % - no detrending, Rectangular window
75 results = [results utp_16]; % Test "conservation of energy":
76 % - white noise from uniform pdf, fixed parameters
77 % - no detrending, Rectangular window
78 results = [results utp_17]; % Test units handling: PSD
79 results = [results utp_18]; % Test units handling: ASD
80 results = [results utp_19]; % Test units handling: PS
81 results = [results utp_20]; % Test units handling: AS
82 results = [results utp_21]; % Test "conservation of energy": Welch window + no detrending
83 results = [results utp_22]; % Test "conservation of energy": Welch window + mean detrending
84 results = [results utp_23]; % Test "conservation of energy": Welch window + linear detrending
85 results = [results utp_24]; % Test "conservation of energy": Hanning window + no detrending
86 results = [results utp_25]; % Test "conservation of energy": Hanning window + mean detrending
87 results = [results utp_26]; % Test "conservation of energy": Hanning window + linear detrending
88 results = [results utp_27]; % Test "conservation of energy": Hamming window + no detrending
89 results = [results utp_28]; % Test "conservation of energy": Hamming window + mean detrending
90 results = [results utp_29]; % Test "conservation of energy": Hamming window + linear detrending
91 results = [results utp_30]; % Test "conservation of energy": BH92 window + no detrending
92 results = [results utp_31]; % Test "conservation of energy": BH92 window + mean detrending
93 results = [results utp_32]; % Test "conservation of energy": BH92 window + linear detrending
94 results = [results utp_33]; % Test "conservation of energy": Kaiser200 window + no detrending
95 results = [results utp_34]; % Test "conservation of energy": Kaiser200 window + mean detrending
96 results = [results utp_35]; % Test "conservation of energy": Kaiser200 window + linear detrending
97 results = [results utp_38]; % Test detrending
98 results = [results utp_39]; % Test detrending
99 results = [results utp_40]; % Test detrending
100 results = [results utp_41]; % Test different windows: Rectangular
101 results = [results utp_42]; % Test different windows: BH92
102 results = [results utp_43]; % Test different windows: Hamming
103 results = [results utp_44]; % Test different windows: Hanning
104 results = [results utp_45]; % Test different windows: Bartlett
105 results = [results utp_46]; % Test different windows: Nuttall3
106 results = [results utp_47]; % Test different windows: Kaiser psll = [random]
107 results = [results utp_48]; % Test different windows: Kaiser psll = [default]
108 results = [results utp_49]; % Test different windows: Nuttall4
109 results = [results utp_50]; % Test different windows: SFT3F
110 results = [results utp_51]; % Test number of averages: requested/obtained
111 results = [results utp_52]; % Test number of averages: correct number
112 results = [results utp_53]; % Test number of averages: syntax
113 results = [results utp_54]; % Test "conservation of energy":
114 results = [results utp_60]; % Test 'basic' call: PSD
115 results = [results utp_61]; % Test 'basic' call: ASD
116 results = [results utp_62]; % Test 'basic' call: PS
117 results = [results utp_63]; % Test 'basic' call: AS
118
119
120 disp('Done.');
121 disp('******************************************************');
122
123 elseif nargin == 1 % Check for UTP functions
124 if strcmp(varargin{1}, 'isutp')
125 results = 1;
126 else
127 results = 0;
128 end
129 else
130 error('### Incorrect inputs')
131 end
132
133
134 %% UTP_01
135
136 % <TestDescription>
137 %
138 % Tests that the getInfo call works for this method.
139 %
140 % </TestDescription>
141 function result = utp_01
142
143
144 % <SyntaxDescription>
145 %
146 % Test that the getInfo call works for no sets, all sets, and each set
147 % individually.
148 %
149 % </SyntaxDescription>
150
151 % <SyntaxCode>
152 try
153 % Call for no sets
154 io(1) = eval([class '.getInfo(''' mthd ''', ''None'')']);
155 % Call for all sets
156 io(2) = eval([class '.getInfo(''' mthd ''')']);
157 % Call for each set
158 for kk=1:numel(io(2).sets)
159 io(kk+2) = eval([class '.getInfo(''' mthd ''', ''' io(2).sets{kk} ''')']);
160 end
161 stest = true;
162 catch err
163 disp(err.message)
164 stest = false;
165 end
166 % </SyntaxCode>
167
168 % <AlgoDescription>
169 %
170 % 1) Check that getInfo call returned an minfo object in all cases.
171 % 2) Check that all plists have the correct parameters.
172 %
173 % </AlgoDescription>
174
175 % <AlgoCode>
176 atest = true;
177 if stest
178 % check we have minfo objects
179 if isa(io, 'minfo')
180
181 %%% SET 'None'
182 if ~isempty(io(1).sets), atest = false; end
183 if ~isempty(io(1).plists), atest = false; end
184 %%% Check all Sets
185 if ~any(strcmpi(io(2).sets, 'Default')), atest = false; end
186 if numel(io(2).plists) ~= numel(io(2).sets), atest = false; end
187 %%%%%%%%%% SET 'Default'
188 if io(3).plists.nparams ~= 9, atest = false; end
189 % Check key
190 if ~io(3).plists.isparam('nfft'), atest = false; end
191 if ~io(3).plists.isparam('win'), atest = false; end
192 if ~io(3).plists.isparam('olap'), atest = false; end
193 if ~io(3).plists.isparam('scale'), atest = false; end
194 if ~io(3).plists.isparam('order'), atest = false; end
195 if ~io(3).plists.isparam('navs'), atest = false; end
196 if ~io(3).plists.isparam('times'), atest = false; end
197 if ~io(3).plists.isparam('split'), atest = false; end
198 if ~io(3).plists.isparam('psll'), atest = false; end
199 % Check default value
200 if ~isequal(io(3).plists.find('nfft'), -1), atest = false; end
201 if ~strcmpi(io(3).plists.find('win'), defaultWinType), atest = false; end
202 if ~isequal(io(3).plists.find('olap'), -1), atest = false; end
203 if ~isequal(io(3).plists.find('scale'), 'PSD'), atest = false; end
204 if ~isequal(io(3).plists.find('order'), 0), atest = false; end
205 if ~isequal(io(3).plists.find('navs'), -1), atest = false; end
206 if ~isEmptyDouble(io(3).plists.find('times')), atest = false; end
207 if ~isEmptyDouble(io(3).plists.find('split')), atest = false; end
208 if ~isequal(io(3).plists.find('psll'), 200), atest = false; end
209 % Check options
210 if ~isequal(io(3).plists.getOptionsForParam('nfft'), {-1}), atest = false; end
211 if ~isequal(io(3).plists.getOptionsForParam('win'), specwin.getTypes), atest = false; end
212 if ~isequal(io(3).plists.getOptionsForParam('olap'), {-1}), atest = false; end
213 if ~isequal(io(3).plists.getOptionsForParam('scale'), {'PSD', 'ASD', 'PS', 'AS'}), atest = false; end
214 if ~isequal(io(3).plists.getOptionsForParam('order'), {-1 0 1 2 3 4 5 6 7 8 9}), atest = false; end
215 if ~isequal(io(3).plists.getOptionsForParam('navs'), {-1}), atest = false; end
216 if ~isequal(io(3).plists.getOptionsForParam('times'), {[]}), atest = false; end
217 if ~isequal(io(3).plists.getOptionsForParam('split'), {[]}), atest = false; end
218 if ~isequal(io(3).plists.getOptionsForParam('psll'), {200}), atest = false; end
219 end
220 else
221 atest = false;
222 end
223 % </AlgoCode>
224
225 % Return a result structure
226 result = utp_prepare_result(atest, stest, dbstack, mfilename);
227 end % END UTP_01
228
229 %% UTP_02
230
231 % <TestDescription>
232 %
233 % Tests that the psd method works with a vector of AOs as input.
234 %
235 % </TestDescription>
236 function result = utp_02
237
238 % <SyntaxDescription>
239 %
240 % Test that the psd method works for a vector of AOs as input.
241 %
242 % </SyntaxDescription>
243
244 % <SyntaxCode>
245 try
246 avec = [at1 at5 at6];
247 % Vector output
248 out = psd(avec);
249 % List output
250 [out1, out2, out3] = psd(avec);
251 stest = true;
252 catch err
253 disp(err.message)
254 stest = false;
255 end
256 % </SyntaxCode>
257
258 % <AlgoDescription>
259 %
260 % 1) Check that the number of elements in 'out' is the same as in the input.
261 % 2) Check that each output object contains the correct values.
262 %
263 % </AlgoDescription>
264
265 % <AlgoCode>
266 atest = true;
267 if stest
268 % Check we have the correct number of outputs
269 if numel(out) ~= numel(avec), atest = false; end
270 % Check we have the correct values in the outputs
271 % Vector output
272 for kk = 1:numel(out)
273 if ~eq(out(kk), psd(avec(kk)), ple1), atest = false; end
274 end
275 % List output
276 if ~eq(out1, psd(avec(1)), ple1), atest = false; end
277 if ~eq(out2, psd(avec(2)), ple1), atest = false; end
278 if ~eq(out3, psd(avec(3)), ple1), atest = false; end
279 else
280 atest = false;
281 end
282 % </AlgoCode>
283
284 % Return a result structure
285 result = utp_prepare_result(atest, stest, dbstack, mfilename);
286 end % END UTP_02
287
288 %% UTP_03
289
290 % <TestDescription>
291 %
292 % Tests that the psd method works with a matrix of AOs as input.
293 %
294 % </TestDescription>
295 function result = utp_03
296
297 % <SyntaxDescription>
298 %
299 % Test that the psd method works for a matrix of AOs as input.
300 %
301 % </SyntaxDescription>
302
303 % <SyntaxCode>
304 try
305 amat = [at1 at5 at6; at5 at6 at1];
306 % Vector output
307 out = psd(amat);
308 % List output
309 [out1, out2, out3, out4, out5, out6] = psd(amat);
310 stest = true;
311 catch err
312 disp(err.message)
313 stest = false;
314 end
315 % </SyntaxCode>
316
317 % <AlgoDescription>
318 %
319 % 1) Check that the number of elements in 'out' is the same as in the input.
320 % 2) Check that each output object contains the correct values.
321 %
322 % </AlgoDescription>
323
324 % <AlgoCode>
325 atest = true;
326 if stest
327 % Check we have the correct number of outputs
328 if numel(out) ~= numel(amat), atest = false; end
329 % Check we have the correct values in the outputs
330 % Vector output
331 for kk = 1:numel(out)
332 if ~eq(out(kk), psd(amat(kk)), ple1), atest = false; end
333 end
334 % List output
335 if ~eq(out1, psd(amat(1)), ple1), atest = false; end
336 if ~eq(out2, psd(amat(2)), ple1), atest = false; end
337 if ~eq(out3, psd(amat(3)), ple1), atest = false; end
338 if ~eq(out4, psd(amat(4)), ple1), atest = false; end
339 if ~eq(out5, psd(amat(5)), ple1), atest = false; end
340 if ~eq(out6, psd(amat(6)), ple1), atest = false; end
341 else
342 atest = false;
343 end
344 % </AlgoCode>
345
346 % Return a result structure
347 result = utp_prepare_result(atest, stest, dbstack, mfilename);
348 end % END UTP_03
349
350 %% UTP_04
351
352 % <TestDescription>
353 %
354 % Tests that the psd method works with a list of AOs as input.
355 %
356 % </TestDescription>
357 function result = utp_04
358
359 % <SyntaxDescription>
360 %
361 % Test that the psd method works for a list of AOs as input.
362 %
363 % </SyntaxDescription>
364
365 % <SyntaxCode>
366 try
367 % Vector output
368 out = psd(at1,at5,at6);
369 % List output
370 [out1, out2, out3] = psd(at1,at5,at6);
371 stest = true;
372 catch err
373 disp(err.message)
374 stest = false;
375 end
376 % </SyntaxCode>
377
378 % <AlgoDescription>
379 %
380 % 1) Check that the number of elements in 'out' is the same as in the input.
381 % 2) Check that each output AO contains the correct data.
382 %
383 % </AlgoDescription>
384
385 % <AlgoCode>
386 atest = true;
387 if stest
388 % Check we have the correct number of outputs
389 if numel(out) ~= 3, atest = false; end
390 % Check we have the correct values in the outputs
391 % Vector output
392 if ~eq(out(1), psd(at1), ple1), atest = false; end
393 if ~eq(out(2), psd(at5), ple1), atest = false; end
394 if ~eq(out(3), psd(at6), ple1), atest = false; end
395 % List output
396 if ~eq(out1, psd(at1), ple1), atest = false; end
397 if ~eq(out2, psd(at5), ple1), atest = false; end
398 if ~eq(out3, psd(at6), ple1), atest = false; end
399
400 else
401 atest = false;
402 end
403 % </AlgoCode>
404
405 % Return a result structure
406 result = utp_prepare_result(atest, stest, dbstack, mfilename);
407 end % END UTP_04
408
409 %% UTP_05
410
411 % <TestDescription>
412 %
413 % Tests that the psd method works with a mix of different shaped AOs as
414 % input.
415 %
416 % </TestDescription>
417 function result = utp_05
418
419 % <SyntaxDescription>
420 %
421 % Test that the psd method works with an input of matrices and vectors
422 % and single AOs.
423 %
424 % </SyntaxDescription>
425
426 % <SyntaxCode>
427 try
428 % Vector output
429 out = psd(at1,[at5 at6],at5,[at5 at1; at6 at1],at6);
430 % List output
431 [out1, out2, out3, out4, out5, out6, out7, out8, out9] = ...
432 psd(at1,[at5 at6],at5,[at5 at1; at6 at1],at6);
433 stest = true;
434 catch err
435 disp(err.message)
436 stest = false;
437 end
438 % </SyntaxCode>
439
440 % <AlgoDescription>
441 %
442 % 1) Check that the number of elements in 'out' is the same as in
443 % input.
444 % 2) Check that each output AO contains the correct data.
445 %
446 % </AlgoDescription>
447
448 % <AlgoCode>
449 atest = true;
450 if stest
451 % Check we have the correct number of outputs
452 if numel(out) ~= 9, atest = false; end
453 % Check we have the correct values in the outputs
454 % Vector output
455 if ~eq(out(1), psd(at1), ple1), atest = false; end
456 if ~eq(out(2), psd(at5), ple1), atest = false; end
457 if ~eq(out(3), psd(at6), ple1), atest = false; end
458 if ~eq(out(4), psd(at5), ple1), atest = false; end
459 if ~eq(out(5), psd(at5), ple1), atest = false; end
460 if ~eq(out(6), psd(at6), ple1), atest = false; end
461 if ~eq(out(7), psd(at1), ple1), atest = false; end
462 if ~eq(out(8), psd(at1), ple1), atest = false; end
463 if ~eq(out(9), psd(at6), ple1), atest = false; end
464 % List output
465 if ~eq(out1, psd(at1), ple1), atest = false; end
466 if ~eq(out2, psd(at5), ple1), atest = false; end
467 if ~eq(out3, psd(at6), ple1), atest = false; end
468 if ~eq(out4, psd(at5), ple1), atest = false; end
469 if ~eq(out5, psd(at5), ple1), atest = false; end
470 if ~eq(out6, psd(at6), ple1), atest = false; end
471 if ~eq(out7, psd(at1), ple1), atest = false; end
472 if ~eq(out8, psd(at1), ple1), atest = false; end
473 if ~eq(out9, psd(at6), ple1), atest = false; end
474
475 else
476 atest = false;
477 end
478 % </AlgoCode>
479
480 % Return a result structure
481 result = utp_prepare_result(atest, stest, dbstack, mfilename);
482 end % END UTP_05
483
484 %% UTP_06
485
486 % <TestDescription>
487 %
488 % Tests that the psd method properly applies history.
489 %
490 % </TestDescription>
491 function result = utp_06
492
493 % <SyntaxDescription>
494 %
495 % Test that the result of applying the psd method can be processed back
496 % to an m-file.
497 %
498 % </SyntaxDescription>
499
500 % <SyntaxCode>
501 try
502 out = psd(at5);
503 mout = rebuild(out);
504 stest = true;
505 catch err
506 disp(err.message)
507 stest = false;
508 end
509 % </SyntaxCode>
510
511 % <AlgoDescription>
512 %
513 % 1) Check that the last entry in the history of 'out' corresponds to
514 % 'psd'.
515 % 2) Check that the re-built object is the same object as 'out'.
516 %
517 % </AlgoDescription>
518
519 % <AlgoCode>
520 atest = true;
521 if stest
522 % Check the last step in the history of 'out'
523 if ~strcmp(out.hist.methodInfo.mname, 'psd'), atest = false; end
524 % Check the re-built object
525 if ~eq(mout, out, ple2), atest = false; end
526 else
527 atest = false;
528 end
529 % </AlgoCode>
530
531 % Return a result structure
532 result = utp_prepare_result(atest, stest, dbstack, mfilename);
533 end % END UTP_06
534
535 %% UTP_07
536
537 % <TestDescription>
538 %
539 % Tests that the psd method agrees with MATLAB's pwelch, within some tolerance,
540 % when configured to use the same parameters.
541 %
542 % </TestDescription>
543 function result = utp_07
544
545 % <SyntaxDescription>
546 %
547 % Test that applying psd works on a single AO.
548 %
549 % </SyntaxDescription>
550
551 % <SyntaxCode>
552 try
553 % Load reference ao
554 a1 = ao('ref_whitenoise_10s_1000Hz.xml');
555 fs = a1.fs;
556 % Compute PSD
557 Nfft = 2*fs;
558 win = specwin('Hanning', Nfft);
559 pl = plist('Nfft', Nfft, 'Win', win.type, 'order', -1);
560 S1 = a1.psd(pl);
561 stest = true;
562 catch err
563 disp(err.message)
564 stest = false;
565 end
566 % </SyntaxCode>
567
568 % <AlgoDescription>
569 %
570 % 1) Check that output agrees with the output of MATLAB's pwelch within tolerance.
571 %
572 % </AlgoDescription>
573
574 % <AlgoCode>
575 atest = true;
576 TOL = 1.8e-15;
577
578 if stest
579 % Load reference psd
580 ref = load('ref_psd_10s_1000Hz.txt');
581 if ~isequal(ref(:,1), S1.x) || any((abs(S1.y - ref(:,2)) ./ S1.y) >= TOL)
582 atest = false;
583 end
584 else
585 atest = false;
586 end
587 % </AlgoCode>
588
589 % Return a result structure
590 result = utp_prepare_result(atest, stest, dbstack, mfilename);
591 end % END UTP_07
592
593 %% UTP_08
594
595 % <TestDescription>
596 %
597 % Check that the psd method pass back the output objects to a list of
598 % output variables or to a single variable.
599 %
600 % </TestDescription>
601 function result = utp_08
602
603 % <SyntaxDescription>
604 %
605 % Call the method with a list of output variables and with a single output
606 % variable. Additionaly check that the rebuild method works on the output.
607 %
608 % </SyntaxDescription>
609
610 try
611 % <SyntaxCode>
612 [o1, o2] = psd(at5, at6);
613 o3 = psd(at5, at6);
614 mout1 = rebuild(o1);
615 mout2 = rebuild(o2);
616 mout3 = rebuild(o3);
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 the output contains the right number of objects
627 % 2) Check that the 'rebuild' method produces the same object as 'out'.
628 %
629 % </AlgoDescription>
630
631 atest = true;
632 if stest
633 % <AlgoCode>
634 % Check the number of outputs
635 if numel(o1) ~=1, atest = false; end
636 if numel(o2) ~=1, atest = false; end
637 if numel(o3) ~=2, atest = false; end
638 % Check the rebuilding of the object
639 if ~eq(o1, mout1, ple2), atest = false; end
640 if ~eq(o2, mout2, ple2), atest = false; end
641 if ~eq(o3, mout3, ple2), atest = false; end
642 % </AlgoCode>
643 else
644 atest = false;
645 end
646
647 % Return a result structure
648 result = utp_prepare_result(atest, stest, dbstack, mfilename);
649 end % END UTP_08
650
651 %% UTP_09
652
653 % <TestDescription>
654 %
655 % Tests that the psd method agrees with MATLAB's pwelch when
656 % configured to use the same parameters.
657 %
658 % </TestDescription>
659 function result = utp_09
660
661 % <SyntaxDescription>
662 %
663 % Test that the applying psd works on a single AO.
664 %
665 % </SyntaxDescription>
666
667 % <SyntaxCode>
668 try
669 % Load reference ao
670 a1 = ao('ref_whitenoise_3600s_1Hz.xml');
671 % Compute PSD
672 Nfft = 1000;
673 win = specwin('BH92', Nfft);
674 pl = plist('Nfft', Nfft, 'Win', win.type, 'order', -1);
675 S1 = a1.psd(pl);
676 stest = true;
677 catch err
678 disp(err.message)
679 stest = false;
680 end
681 % </SyntaxCode>
682
683 % <AlgoDescription>
684 %
685 % 1) Check that output agrees with the output of MATLAB's pwelch.
686 %
687 % </AlgoDescription>
688
689 % <AlgoCode>
690 atest = true;
691 TOL = 1e-15;
692
693 if stest
694 % Load reference psd
695 ref = load('ref_psd_3600s_1Hz.txt');
696 if ~isequal(ref(:,1), S1.x) || any((abs(S1.y - ref(:,2)) ./ S1.y) >= TOL)
697 atest = false;
698 end
699 else
700 atest = false;
701 end
702 % </AlgoCode>
703
704 % Return a result structure
705 result = utp_prepare_result(atest, stest, dbstack, mfilename);
706 end % END UTP_09
707
708 %% UTP_10
709
710 % <TestDescription>
711 %
712 % Tests that the psd method agrees with MATLAB's pwelch when
713 % configured to use the same parameters.
714 %
715 % </TestDescription>
716 function result = utp_10
717
718 % <SyntaxDescription>
719 %
720 % Test that the applying psd works on a single AO.
721 %
722 % </SyntaxDescription>
723
724 % <SyntaxCode>
725 try
726 % Load reference ao
727 a1 = ao('ref_whitenoise_100000s_100mHz.xml');
728 % Compute PSD
729 Nfft = a1.fs * a1.data.nsecs;
730 psll = 100;
731 win = specwin('Kaiser', Nfft, psll);
732 pl = plist('Nfft', Nfft, 'Win', win.type, 'psll', psll, 'order', -1);
733 S1 = a1.psd(pl);
734 stest = true;
735 catch err
736 disp(err.message)
737 stest = false;
738 end
739 % </SyntaxCode>
740
741 % <AlgoDescription>
742 %
743 % 1) Check that output agrees with the output of MATLAB's pwelch.
744 %
745 % </AlgoDescription>
746
747 % <AlgoCode>
748 atest = true;
749 TOL = 1e-13;
750
751 if stest
752 % Load reference psd
753 ref = load('ref_psd_100000s_100mHz.txt');
754 if ~isequal(ref(:,1), S1.x) || any((abs(S1.y - ref(:,2)) ./ S1.y) >= TOL)
755 atest = false;
756 end
757 else
758 atest = false;
759 end
760 % </AlgoCode>
761
762 % Return a result structure
763 result = utp_prepare_result(atest, stest, dbstack, mfilename);
764 end % END UTP_10
765
766
767
768 %% UTP_12
769
770 % <TestDescription>
771 %
772 % Tests "conservation of energy":
773 % 1) white noise produced from uniform pdf, with a given mean value and
774 % sigma (distribution's 1st and 2nd order)
775 % 2) evaluate the sample mean value m and standard deviation s
776 % 3) psd of the white noise
777 % 4) compares the sqrt(sum(S * df)) with the standard deviation s
778 %
779
780 % </TestDescription>
781 function result = utp_12
782
783 % <SyntaxDescription>
784 %
785 % 1) Prepare the test tsdata:
786 % white noise from uniform distribution + offset
787 % 2) Calculate the statistical parameters
788 % 3) Estimate the psd
789 %
790 % </SyntaxDescription>
791
792 % <SyntaxCode>
793 try
794
795 % White noise
796 % Parameters are crucial for the estimate to be correct!
797
798 noise_type = 'Uniform';
799 win_type = 'Rectangular';
800 olap = 0;
801 detrend_order = 0;
802
803 scale = 'PSD';
804
805 [sigma_sample, sigma_psd] = prepare_analyze_noise_energy(win_type, noise_type, scale, ...
806 plist('detrend_order', detrend_order, 'olap', olap));
807
808 stest = true;
809
810 catch err
811 disp(err.message)
812 stest = false;
813 end
814
815 % </SyntaxCode>
816
817 % <AlgoDescription>
818 %
819 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
820 %
821 % </AlgoDescription>
822
823 % <AlgoCode>
824 atest = true;
825 TOL = 1e-12;
826
827 if stest
828 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL
829 atest = false;
830 end
831 else
832 atest = false;
833 end
834 % </AlgoCode>
835
836 % Return a result structure
837 result = utp_prepare_result(atest, stest, dbstack, mfilename);
838 end % END UTP_12
839
840 %% UTP_13
841
842 % <TestDescription>
843 %
844 % Tests "conservation of energy" with fixed parameters:
845 % 1) white noise produced from normal pdf, with a given mean value and
846 % sigma (distribution's 1st and 2nd order)
847 % 2) evaluate the sample mean value m and standard deviation s
848 % 3) psd of the white noise
849 % 4) compares the sqrt(sum(S * df)) with the standard deviation s
850 %
851
852 % </TestDescription>
853 function result = utp_13
854
855 % <SyntaxDescription>
856 %
857 % 1) Prepare the test tsdata:
858 % white noise from normal distribution + offset
859 % 2) Calculate the statistical parameters
860 % 3) Estimate the psd
861 %
862 % </SyntaxDescription>
863
864 % <SyntaxCode>
865 try
866
867 % White noise
868 % Parameters are crucial for the estimate to be correct!
869
870 noise_type = 'Normal';
871 win_type = 'Rectangular';
872 olap = 0;
873 detrend_order = 0;
874
875 scale = 'PSD';
876
877 % Build with fixed parameters
878 fs = 10;
879 nsecs = 3600;
880 sigma_distr = 1;
881 mu_distr = 0;
882 [sigma_sample, sigma_psd] = prepare_analyze_noise_energy(win_type, noise_type, scale, ...
883 plist('detrend_order', detrend_order, 'olap', olap, ...
884 'fs', fs, 'nsecs', nsecs, 'sigma_distr', sigma_distr, 'mu_distr', mu_distr));
885
886 stest = true;
887
888 catch err
889 disp(err.message)
890 stest = false;
891 end
892 % </SyntaxCode>
893
894 % <AlgoDescription>
895 %
896 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
897 %
898 % </AlgoDescription>
899
900 % <AlgoCode>
901 atest = true;
902 TOL = 1e-12;
903
904 if stest
905 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL
906 atest = false;
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 "conservation of energy" with fixed parameters:
922 % 1) white noise produced from uniform pdf, with a given mean value and
923 % sigma (distribution's 1st and 2nd order)
924 % 2) evaluate the sample mean value m and standard deviation s
925 % 3) psd of the white noise
926 % 4) compares the sqrt(sum(S * df)) with the standard deviation s
927 %
928
929 % </TestDescription>
930 function result = utp_14
931
932 % <SyntaxDescription>
933 %
934 % 1) Prepare the test tsdata:
935 % white noise from uniform distribution + offset
936 % 2) Calculate the statistical parameters
937 % 3) Estimate the psd
938 %
939 % </SyntaxDescription>
940
941 % <SyntaxCode>
942 try
943
944 % White noise
945 % Parameters are crucial for the estimate to be correct!
946
947 noise_type = 'Uniform';
948 win_type = 'Rectangular';
949 olap = 0;
950 detrend_order = 0;
951
952 scale = 'PSD';
953
954 % Build with fixed parameters
955 fs = 1;
956 nsecs = 86400;
957 sigma_distr = 1e-9;
958 mu_distr = 3e-7;
959 [sigma_sample, sigma_psd] = prepare_analyze_noise_energy(win_type, noise_type, scale, ...
960 plist('detrend_order', detrend_order, 'olap', olap, ...
961 'fs', fs, 'nsecs', nsecs, 'sigma_distr', sigma_distr, 'mu_distr', mu_distr));
962
963 stest = true;
964
965 catch err
966 disp(err.message)
967 stest = false;
968 end
969 % </SyntaxCode>
970
971 % <AlgoDescription>
972 %
973 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
974 %
975 % </AlgoDescription>
976
977 % <AlgoCode>
978 atest = true;
979 TOL = 1e-12;
980
981 if stest
982 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL
983 atest = false;
984 end
985 else
986 atest = false;
987 end
988 % </AlgoCode>
989
990 % Return a result structure
991 result = utp_prepare_result(atest, stest, dbstack, mfilename);
992 end % END UTP_14
993
994 %% UTP_15
995
996 % <TestDescription>
997 %
998 % Tests "conservation of energy" with fixed parameters:
999 % 1) white noise produced from normal pdf, with a given mean value and
1000 % sigma (distribution's 1st and 2nd order)
1001 % 2) evaluate the sample mean value m and standard deviation s
1002 % 3) psd of the white noise
1003 % 4) compares the sqrt(sum(S * df)) with the standard deviation s
1004 %
1005
1006 % </TestDescription>
1007 function result = utp_15
1008
1009 % <SyntaxDescription>
1010 %
1011 % 1) Prepare the test tsdata:
1012 % white noise from normal distribution + offset
1013 % 2) Calculate the statistical parameters
1014 % 3) Estimate the psd
1015 %
1016 % </SyntaxDescription>
1017
1018 % <SyntaxCode>
1019 try
1020
1021 % White noise
1022 % Parameters are crucial for the estimate to be correct!
1023
1024 noise_type = 'Normal';
1025 win_type = 'Rectangular';
1026 olap = 0;
1027 detrend_order = 0;
1028
1029 scale = 'PSD';
1030
1031 % Build with fixed parameters
1032 fs = 100;
1033 nsecs = 100;
1034 sigma_distr = 1e-12;
1035 mu_distr = -5e-12;
1036 [sigma_sample, sigma_psd] = prepare_analyze_noise_energy(win_type, noise_type, scale, ...
1037 plist('detrend_order', detrend_order, 'olap', olap, ...
1038 'fs', fs, 'nsecs', nsecs, 'sigma_distr', sigma_distr, 'mu_distr', mu_distr));
1039
1040 stest = true;
1041
1042 catch err
1043 disp(err.message)
1044 stest = false;
1045 end
1046 % </SyntaxCode>
1047
1048 % <AlgoDescription>
1049 %
1050 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
1051 %
1052 % </AlgoDescription>
1053
1054 % <AlgoCode>
1055 atest = true;
1056 TOL = 1e-12;
1057
1058 if stest
1059 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL
1060 atest = false;
1061 end
1062 else
1063 atest = false;
1064 end
1065 % </AlgoCode>
1066
1067 % Return a result structure
1068 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1069 end % END UTP_15
1070
1071 %% UTP_16
1072
1073 % <TestDescription>
1074 %
1075 % Tests "conservation of energy" with fixed parameters:
1076 % 1) white noise produced from uniform pdf, with a given mean value and
1077 % sigma (distribution's 1st and 2nd order)
1078 % 2) evaluate the sample mean value m and standard deviation s
1079 % 3) psd of the white noise
1080 % 4) compares the sqrt(sum(S * df)) with the standard deviation s
1081 %
1082
1083 % </TestDescription>
1084 function result = utp_16
1085
1086 % <SyntaxDescription>
1087 %
1088 % 1) Prepare the test tsdata:
1089 % white noise from uniform distribution + offset
1090 % 2) Calculate the statistical parameters
1091 % 3) Estimate the psd
1092 %
1093 % </SyntaxDescription>
1094
1095 % <SyntaxCode>
1096 try
1097
1098 % White noise
1099 % Parameters are crucial for the estimate to be correct!
1100
1101 noise_type = 'Uniform';
1102 win_type = 'Rectangular';
1103 olap = 0;
1104 detrend_order = 0;
1105
1106 scale = 'PSD';
1107
1108 % Build with fixed parameters
1109 fs = 0.01;
1110 nsecs = 3*86400;
1111 sigma_distr = 1e-15;
1112 mu_distr = -7.72e-13;
1113 [sigma_sample, sigma_psd] = prepare_analyze_noise_energy(win_type, noise_type, scale, ...
1114 plist('detrend_order', detrend_order, 'olap', olap, ...
1115 'fs', fs, 'nsecs', nsecs, 'sigma_distr', sigma_distr, 'mu_distr', mu_distr));
1116
1117 stest = true;
1118
1119 catch err
1120 disp(err.message)
1121 stest = false;
1122 end
1123 % </SyntaxCode>
1124
1125 % <AlgoDescription>
1126 %
1127 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
1128 %
1129 % </AlgoDescription>
1130
1131 % <AlgoCode>
1132 atest = true;
1133 TOL = 1e-12;
1134
1135 if stest
1136 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL
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_16
1147
1148 %% UTP_17
1149
1150 % <TestDescription>
1151 %
1152 % Tests handling of units:
1153 % 1) white noise produced from normal pdf, with a given mean value and
1154 % sigma (distribution's 1st and 2nd orders)
1155 % 2) PSD of the white noise
1156 % 3) compares the units of the input and output
1157 %
1158
1159 % </TestDescription>
1160 function result = utp_17
1161
1162 % <SyntaxDescription>
1163 %
1164 % 1) Prepare the test tsdata:
1165 % white noise from normal distribution + offset
1166 % 2) Assign a random unit
1167 % 3) PSD of the white noise
1168 %
1169 % </SyntaxDescription>
1170
1171 % <SyntaxCode>
1172 try
1173
1174 % White noise
1175 noise_type = 'Normal';
1176 win_type = 'BH92';
1177 reduce_pts = 10;
1178 scale = 'PSD';
1179 olap = specwin(win_type).rov;
1180
1181 [a, S, S1, S2] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
1182 plist('olap', olap, 'reduce_pts', reduce_pts));
1183 % S1 and S2 are empty in this case
1184
1185 stest = true;
1186
1187 catch err
1188 disp(err.message)
1189 stest = false;
1190 end
1191 % </SyntaxCode>
1192
1193 % <AlgoDescription>
1194 %
1195 % 1) Check that (calculated PSD yunits) equals (input units)^2 / Hz
1196 % 2) Check that (calculated PSD xunits) equals Hz
1197 %
1198 % </AlgoDescription>
1199
1200 % <AlgoCode>
1201 atest = true;
1202
1203 if stest
1204 if ne(S.yunits, (a.yunits).^2 * unit('Hz^-1')) || ne(S.xunits, unit('Hz'))
1205 atest = false;
1206 end
1207 else
1208 atest = false;
1209 end
1210 % </AlgoCode>
1211
1212 % Return a result structure
1213 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1214 end % END UTP_17
1215
1216 %% UTP_18
1217
1218 % <TestDescription>
1219 %
1220 % Tests handling of units:
1221 % 1) white noise produced from uniform pdf, with a given mean value and
1222 % sigma (distribution's 1st and 2nd orders)
1223 % 2) ASD of the white noise
1224 % 3) compares the units of the input and output
1225 %
1226
1227 % </TestDescription>
1228 function result = utp_18
1229
1230 % <SyntaxDescription>
1231 %
1232 % 1) Prepare the test tsdata:
1233 % white noise from uniform distribution + offset
1234 % 2) Assign a random unit
1235 % 3) ASD of the white noise
1236 %
1237 % </SyntaxDescription>
1238
1239 % <SyntaxCode>
1240 try
1241
1242 % White noise
1243 noise_type = 'Uniform';
1244 win_type = 'Hamming';
1245 reduce_pts = 4;
1246 scale = 'ASD';
1247 olap = specwin(win_type).rov;
1248
1249 [a, S, S1, S2] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
1250 plist('olap', olap, 'reduce_pts', reduce_pts));
1251 % S1 and S2 are empty in this case
1252
1253 stest = true;
1254
1255 catch err
1256 disp(err.message)
1257 stest = false;
1258 end
1259 % </SyntaxCode>
1260
1261 % <AlgoDescription>
1262 %
1263 % 1) Check that (calculated ASD yunits) equals (input units) / Hz^(1/2)
1264 % 2) Check that (calculated ASD xunits) equals Hz
1265 %
1266 % </AlgoDescription>
1267
1268 % <AlgoCode>
1269 atest = true;
1270
1271 if stest
1272 if ne(S.yunits, (a.yunits) * unit('Hz^-1/2')) || ne(S.xunits, unit('Hz'))
1273 atest = false;
1274 end
1275 else
1276 atest = false;
1277 end
1278 % </AlgoCode>
1279
1280 % Return a result structure
1281 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1282 end % END UTP_18
1283
1284 %% UTP_19
1285
1286 % <TestDescription>
1287 %
1288 % Tests handling of units:
1289 % 1) white noise produced from normal pdf, with a given mean value and
1290 % sigma (distribution's 1st and 2nd orders)
1291 % 2) PS of the white noise
1292 % 3) compares the units of the input and output
1293 %
1294
1295 % </TestDescription>
1296 function result = utp_19
1297
1298 % <SyntaxDescription>
1299 %
1300 % 1) Prepare the test tsdata:
1301 % white noise from normal distribution + offset
1302 % 2) Assign a random unit
1303 % 3) PS of the white noise
1304 %
1305 % </SyntaxDescription>
1306
1307 % <SyntaxCode>
1308 try
1309
1310 % White noise
1311 noise_type = 'Normal';
1312 win_type = 'Hanning';
1313 reduce_pts = 5;
1314 scale = 'PS';
1315 olap = specwin(win_type).rov;
1316
1317 [a, S, S1, S2] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
1318 plist('olap', olap, 'reduce_pts', reduce_pts));
1319 % S1 and S2 are empty in this case
1320
1321 stest = true;
1322
1323 catch err
1324 disp(err.message)
1325 stest = false;
1326 end
1327 % </SyntaxCode>
1328
1329 % <AlgoDescription>
1330 %
1331 % 1) Check that (calculated PS yunits) equals (input units)^2
1332 % 2) Check that (calculated PS xunits) equals Hz
1333 %
1334 % </AlgoDescription>
1335
1336 % <AlgoCode>
1337 atest = true;
1338
1339 if stest
1340 if ne(S.yunits, (a.yunits).^2) || ne(S.xunits, unit('Hz'))
1341 atest = false;
1342 end
1343 else
1344 atest = false;
1345 end
1346 % </AlgoCode>
1347
1348 % Return a result structure
1349 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1350 end % END UTP_19
1351
1352 %% UTP_20
1353
1354 % <TestDescription>
1355 %
1356 % Tests handling of units:
1357 % 1) white noise produced from uniform distribution, with a given mean value and
1358 % sigma (distribution's 1st and 2nd orders)
1359 % 2) AS of the white noise
1360 % 3) compares the units of the input and output
1361 %
1362
1363 % </TestDescription>
1364 function result = utp_20
1365
1366 % <SyntaxDescription>
1367 %
1368 % 1) Prepare the test tsdata:
1369 % white noise from uniform distribution + offset
1370 % 2) Assign a random unit
1371 % 3) AS of the white noise
1372 %
1373 % </SyntaxDescription>
1374
1375 % <SyntaxCode>
1376 try
1377
1378 % White noise
1379 noise_type = 'Uniform';
1380 win_type = 'Rectangular';
1381 reduce_pts = 10;
1382 scale = 'AS';
1383 olap = specwin(win_type).rov;
1384
1385 [a, S, S1, S2] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
1386 plist('olap', olap, 'reduce_pts', reduce_pts));
1387 % S1 and S2 are empty in this case
1388
1389 stest = true;
1390
1391 catch err
1392 disp(err.message)
1393 stest = false;
1394 end
1395 % </SyntaxCode>
1396
1397 % <AlgoDescription>
1398 %
1399 % 1) Check that (calculated AS yunits) equals (input units)
1400 % 2) Check that (calculated AS xunits) equals Hz
1401 %
1402 % </AlgoDescription>
1403
1404 % <AlgoCode>
1405 atest = true;
1406
1407 if stest
1408 if ne(S.yunits, a.yunits) || ne(S.xunits, unit('Hz'))
1409 atest = false;
1410 end
1411 else
1412 atest = false;
1413 end
1414 % </AlgoCode>
1415
1416 % Return a result structure
1417 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1418 end % END UTP_20
1419
1420 %% UTP_21
1421
1422 % <TestDescription>
1423 %
1424 % Tests "conservation of energy":
1425 % 1) white noise produced from normal pdf, with:
1426 % a given mean value and sigma (distribution's 1st and 2nd order)
1427 % 2) Calculate the expected level of noise from sample statistics
1428 % 3) Calculates the PSD of the individual parts, with:
1429 % Welch window and no detrend
1430 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1431 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
1432 % 6) Evaluate the expected value, estimated from the std of the
1433 % individual segments
1434 % 7) Compares with the mean performed on the individual segments
1435 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1436 % 9) Compares the grand-mean with the estimated white noise level
1437
1438 % </TestDescription>
1439 function result = utp_21
1440
1441 % <SyntaxDescription>
1442 %
1443 % 1) Split the reference data into different segments
1444 % 2) Calculates the PSD of the individual parts
1445 %
1446 % </SyntaxDescription>
1447
1448 % <SyntaxCode>
1449 try
1450
1451 % Calculate the expected level of noise from sample statistics
1452
1453 % Evaluate the PSDs on the parts, employing:
1454 % - Welch window
1455 % - no detrend
1456 win_type = 'Welch';
1457 detrend_order = -1;
1458
1459 [S_ap, S_sigma_total, S_expected_total] = ...
1460 test_psd_energy_prepare_data(win_type, detrend_order, plist);
1461
1462 stest = true;
1463
1464 catch err
1465 disp(err.message)
1466 stest = false;
1467 end
1468 % </SyntaxCode>
1469
1470 % <AlgoDescription>
1471 %
1472 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1473 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
1474 % 3) Evaluate the expected value, estimated from the std of the
1475 % individual segments
1476 % 4) Compares with the mean performed on the individual segments
1477 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1478 % 6) Compares the grand-mean with the estimated white noise level
1479
1480 % </AlgoDescription>
1481
1482 % <AlgoCode>
1483
1484 if stest
1485 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
1486 else
1487 atest = false;
1488 end
1489
1490 % </AlgoCode>
1491
1492 % Return a result structure
1493 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1494 end % END UTP_21
1495
1496
1497 %% UTP_22
1498
1499 % <TestDescription>
1500 %
1501 % Tests "conservation of energy":
1502 % 1) white noise produced from normal pdf, with:
1503 % a given mean value and sigma (distribution's 1st and 2nd order)
1504 % 2) Calculate the expected level of noise from sample statistics
1505 % 3) Calculates the PSD of the individual parts, with:
1506 % Welch window and mean detrend
1507 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1508 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
1509 % 6) Evaluate the expected value, estimated from the std of the
1510 % individual segments
1511 % 7) Compares with the mean performed on the individual segments
1512 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1513 % 9) Compares the grand-mean with the estimated white noise level
1514
1515 % </TestDescription>
1516 function result = utp_22
1517
1518 % <SyntaxDescription>
1519 %
1520 % 1) Split the reference data into different segments
1521 % 2) Calculates the PSD of the individual parts
1522 %
1523 % </SyntaxDescription>
1524
1525 % <SyntaxCode>
1526 try
1527
1528 % Calculate the expected level of noise from sample statistics
1529
1530 % Evaluate the PSDs on the parts, employing:
1531 % - Welch window
1532 % - mean detrend
1533 win_type = 'Welch';
1534 detrend_order = 0;
1535
1536 [S_ap, S_sigma_total, S_expected_total] = ...
1537 test_psd_energy_prepare_data(win_type, detrend_order, plist);
1538
1539 stest = true;
1540
1541 catch err
1542 disp(err.message)
1543 stest = false;
1544 end
1545 % </SyntaxCode>
1546
1547 % <AlgoDescription>
1548 %
1549 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1550 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
1551 % 3) Evaluate the expected value, estimated from the std of the
1552 % individual segments
1553 % 4) Compares with the mean performed on the individual segments
1554 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1555 % 6) Compares the grand-mean with the estimated white noise level
1556
1557 % </AlgoDescription>
1558
1559 % <AlgoCode>
1560
1561 if stest
1562 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
1563 else
1564 atest = false;
1565 end
1566
1567 % </AlgoCode>
1568
1569 % Return a result structure
1570 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1571 end % END UTP_22
1572
1573 %% UTP_23
1574
1575 % <TestDescription>
1576 %
1577 % Tests "conservation of energy":
1578 % 1) white noise produced from normal pdf, with:
1579 % a given mean value and sigma (distribution's 1st and 2nd order)
1580 % 2) Calculate the expected level of noise from sample statistics
1581 % 3) Calculates the PSD of the individual parts, with:
1582 % Welch window and linear detrend
1583 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1584 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
1585 % 6) Evaluate the expected value, estimated from the std of the
1586 % individual segments
1587 % 7) Compares with the mean performed on the individual segments
1588 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1589 % 9) Compares the grand-mean with the estimated white noise level
1590
1591 % </TestDescription>
1592 function result = utp_23
1593
1594 % <SyntaxDescription>
1595 %
1596 % 1) Split the reference data into different segments
1597 % 2) Calculates the PSD of the individual parts
1598 %
1599 % </SyntaxDescription>
1600
1601 % <SyntaxCode>
1602 try
1603
1604 % Calculate the expected level of noise from sample statistics
1605
1606 % Evaluate the PSDs on the parts, employing:
1607 % - Welch window
1608 % - linear detrend
1609 win_type = 'Welch';
1610 detrend_order = 1;
1611
1612 [S_ap, S_sigma_total, S_expected_total] = ...
1613 test_psd_energy_prepare_data(win_type, detrend_order, plist);
1614
1615 stest = true;
1616
1617 catch err
1618 disp(err.message)
1619 stest = false;
1620 end
1621 % </SyntaxCode>
1622
1623 % <AlgoDescription>
1624 %
1625 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1626 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
1627 % 3) Evaluate the expected value, estimated from the std of the
1628 % individual segments
1629 % 4) Compares with the mean performed on the individual segments
1630 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1631 % 6) Compares the grand-mean with the estimated white noise level
1632
1633 % </AlgoDescription>
1634
1635 % <AlgoCode>
1636
1637 if stest
1638 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
1639 else
1640 atest = false;
1641 end
1642
1643 % </AlgoCode>
1644
1645 % Return a result structure
1646 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1647 end % END UTP_23
1648
1649 %% UTP_24
1650
1651 % <TestDescription>
1652 %
1653 % Tests "conservation of energy":
1654 % 1) white noise produced from normal pdf, with:
1655 % a given mean value and sigma (distribution's 1st and 2nd order)
1656 % 2) Calculate the expected level of noise from sample statistics
1657 % 3) Calculates the PSD of the individual parts, with:
1658 % Hanning window and no detrend
1659 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1660 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
1661 % 6) Evaluate the expected value, estimated from the std of the
1662 % individual segments
1663 % 7) Compares with the mean performed on the individual segments
1664 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1665 % 9) Compares the grand-mean with the estimated white noise level
1666
1667 % </TestDescription>
1668 function result = utp_24
1669
1670 % <SyntaxDescription>
1671 %
1672 % 1) Split the reference data into different segments
1673 % 2) Calculates the PSD of the individual parts
1674 %
1675 % </SyntaxDescription>
1676
1677 % <SyntaxCode>
1678 try
1679
1680 % Calculate the expected level of noise from sample statistics
1681
1682 % Evaluate the PSDs on the parts, employing:
1683 % - Hanning window
1684 % - no detrend
1685 win_type = 'Hanning';
1686 detrend_order = -1;
1687
1688 [S_ap, S_sigma_total, S_expected_total] = ...
1689 test_psd_energy_prepare_data(win_type, detrend_order, plist);
1690
1691 stest = true;
1692
1693 catch err
1694 disp(err.message)
1695 stest = false;
1696 end
1697 % </SyntaxCode>
1698
1699 % <AlgoDescription>
1700 %
1701 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1702 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
1703 % 3) Evaluate the expected value, estimated from the std of the
1704 % individual segments
1705 % 4) Compares with the mean performed on the individual segments
1706 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1707 % 6) Compares the grand-mean with the estimated white noise level
1708
1709 % </AlgoDescription>
1710
1711 % <AlgoCode>
1712
1713 if stest
1714 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
1715 else
1716 atest = false;
1717 end
1718
1719 % </AlgoCode>
1720
1721 % Return a result structure
1722 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1723 end % END UTP_24
1724
1725
1726 %% UTP_25
1727
1728 % <TestDescription>
1729 %
1730 % Tests "conservation of energy":
1731 % 1) white noise produced from normal pdf, with:
1732 % a given mean value and sigma (distribution's 1st and 2nd order)
1733 % 2) Calculate the expected level of noise from sample statistics
1734 % 3) Calculates the PSD of the individual parts, with:
1735 % Hanning window and mean detrend
1736 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1737 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
1738 % 6) Evaluate the expected value, estimated from the std of the
1739 % individual segments
1740 % 7) Compares with the mean performed on the individual segments
1741 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1742 % 9) Compares the grand-mean with the estimated white noise level
1743
1744 % </TestDescription>
1745 function result = utp_25
1746
1747 % <SyntaxDescription>
1748 %
1749 % 1) Split the reference data into different segments
1750 % 2) Calculates the PSD of the individual parts
1751 %
1752 % </SyntaxDescription>
1753
1754 % <SyntaxCode>
1755 try
1756
1757 % Calculate the expected level of noise from sample statistics
1758
1759 % Evaluate the PSDs on the parts, employing:
1760 % - Hanning window
1761 % - mean detrend
1762 win_type = 'Hanning';
1763 detrend_order = 0;
1764
1765 [S_ap, S_sigma_total, S_expected_total] = ...
1766 test_psd_energy_prepare_data(win_type, detrend_order, plist);
1767
1768 stest = true;
1769
1770 catch err
1771 disp(err.message)
1772 stest = false;
1773 end
1774 % </SyntaxCode>
1775
1776 % <AlgoDescription>
1777 %
1778 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1779 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
1780 % 3) Evaluate the expected value, estimated from the std of the
1781 % individual segments
1782 % 4) Compares with the mean performed on the individual segments
1783 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1784 % 6) Compares the grand-mean with the estimated white noise level
1785
1786 % </AlgoDescription>
1787
1788 % <AlgoCode>
1789
1790 if stest
1791 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
1792 else
1793 atest = false;
1794 end
1795
1796 % </AlgoCode>
1797
1798 % Return a result structure
1799 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1800 end % END UTP_25
1801
1802 %% UTP_26
1803
1804 % <TestDescription>
1805 %
1806 % Tests "conservation of energy":
1807 % 1) white noise produced from normal pdf, with:
1808 % a given mean value and sigma (distribution's 1st and 2nd order)
1809 % 2) Calculate the expected level of noise from sample statistics
1810 % 3) Calculates the PSD of the individual parts, with:
1811 % Hanning window and linear detrend
1812 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1813 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
1814 % 6) Evaluate the expected value, estimated from the std of the
1815 % individual segments
1816 % 7) Compares with the mean performed on the individual segments
1817 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1818 % 9) Compares the grand-mean with the estimated white noise level
1819
1820 % </TestDescription>
1821 function result = utp_26
1822
1823 % <SyntaxDescription>
1824 %
1825 % 1) Split the reference data into different segments
1826 % 2) Calculates the PSD of the individual parts
1827 %
1828 % </SyntaxDescription>
1829
1830 % <SyntaxCode>
1831 try
1832
1833 % Calculate the expected level of noise from sample statistics
1834
1835 % Evaluate the PSDs on the parts, employing:
1836 % - Hanning window
1837 % - linear detrend
1838 win_type = 'Hanning';
1839 detrend_order = 1;
1840
1841 [S_ap, S_sigma_total, S_expected_total] = ...
1842 test_psd_energy_prepare_data(win_type, detrend_order, plist);
1843
1844 stest = true;
1845
1846 catch err
1847 disp(err.message)
1848 stest = false;
1849 end
1850 % </SyntaxCode>
1851
1852 % <AlgoDescription>
1853 %
1854 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1855 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
1856 % 3) Evaluate the expected value, estimated from the std of the
1857 % individual segments
1858 % 4) Compares with the mean performed on the individual segments
1859 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1860 % 6) Compares the grand-mean with the estimated white noise level
1861
1862 % </AlgoDescription>
1863
1864 % <AlgoCode>
1865
1866 if stest
1867 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
1868 else
1869 atest = false;
1870 end
1871
1872 % </AlgoCode>
1873
1874 % Return a result structure
1875 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1876 end % END UTP_26
1877
1878 %% UTP_27
1879
1880 % <TestDescription>
1881 %
1882 % Tests "conservation of energy":
1883 % 1) white noise produced from normal pdf, with:
1884 % a given mean value and sigma (distribution's 1st and 2nd order)
1885 % 2) Calculate the expected level of noise from sample statistics
1886 % 3) Calculates the PSD of the individual parts, with:
1887 % Hamming window and no detrend
1888 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1889 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
1890 % 6) Evaluate the expected value, estimated from the std of the
1891 % individual segments
1892 % 7) Compares with the mean performed on the individual segments
1893 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1894 % 9) Compares the grand-mean with the estimated white noise level
1895
1896 % </TestDescription>
1897 function result = utp_27
1898
1899 % <SyntaxDescription>
1900 %
1901 % 1) Split the reference data into different segments
1902 % 2) Calculates the PSD of the individual parts
1903 %
1904 % </SyntaxDescription>
1905
1906 % <SyntaxCode>
1907 try
1908
1909 % Calculate the expected level of noise from sample statistics
1910
1911 % Evaluate the PSDs on the parts, employing:
1912 % - Hamming window
1913 % - no detrend
1914 win_type = 'Hamming';
1915 detrend_order = -1;
1916
1917 [S_ap, S_sigma_total, S_expected_total] = ...
1918 test_psd_energy_prepare_data(win_type, detrend_order, plist);
1919
1920 stest = true;
1921
1922 catch err
1923 disp(err.message)
1924 stest = false;
1925 end
1926 % </SyntaxCode>
1927
1928 % <AlgoDescription>
1929 %
1930 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1931 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
1932 % 3) Evaluate the expected value, estimated from the std of the
1933 % individual segments
1934 % 4) Compares with the mean performed on the individual segments
1935 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1936 % 6) Compares the grand-mean with the estimated white noise level
1937
1938 % </AlgoDescription>
1939
1940 % <AlgoCode>
1941
1942 if stest
1943 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
1944 else
1945 atest = false;
1946 end
1947
1948 % </AlgoCode>
1949
1950 % Return a result structure
1951 result = utp_prepare_result(atest, stest, dbstack, mfilename);
1952 end % END UTP_27
1953
1954
1955 %% UTP_28
1956
1957 % <TestDescription>
1958 %
1959 % Tests "conservation of energy":
1960 % 1) white noise produced from normal pdf, with:
1961 % a given mean value and sigma (distribution's 1st and 2nd order)
1962 % 2) Calculate the expected level of noise from sample statistics
1963 % 3) Calculates the PSD of the individual parts, with:
1964 % Hamming window and mean detrend
1965 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
1966 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
1967 % 6) Evaluate the expected value, estimated from the std of the
1968 % individual segments
1969 % 7) Compares with the mean performed on the individual segments
1970 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
1971 % 9) Compares the grand-mean with the estimated white noise level
1972
1973 % </TestDescription>
1974 function result = utp_28
1975
1976 % <SyntaxDescription>
1977 %
1978 % 1) Split the reference data into different segments
1979 % 2) Calculates the PSD of the individual parts
1980 %
1981 % </SyntaxDescription>
1982
1983 % <SyntaxCode>
1984 try
1985
1986 % Calculate the expected level of noise from sample statistics
1987
1988 % Evaluate the PSDs on the parts, employing:
1989 % - Hamming window
1990 % - mean detrend
1991 win_type = 'Hanning';
1992 detrend_order = 0;
1993
1994 [S_ap, S_sigma_total, S_expected_total] = ...
1995 test_psd_energy_prepare_data(win_type, detrend_order, plist);
1996
1997 stest = true;
1998
1999 catch err
2000 disp(err.message)
2001 stest = false;
2002 end
2003 % </SyntaxCode>
2004
2005 % <AlgoDescription>
2006 %
2007 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2008 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
2009 % 3) Evaluate the expected value, estimated from the std of the
2010 % individual segments
2011 % 4) Compares with the mean performed on the individual segments
2012 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2013 % 6) Compares the grand-mean with the estimated white noise level
2014
2015 % </AlgoDescription>
2016
2017 % <AlgoCode>
2018
2019 if stest
2020 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
2021 else
2022 atest = false;
2023 end
2024
2025 % </AlgoCode>
2026
2027 % Return a result structure
2028 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2029 end % END UTP_28
2030
2031 %% UTP_29
2032
2033 % <TestDescription>
2034 %
2035 % Tests "conservation of energy":
2036 % 1) white noise produced from normal pdf, with:
2037 % a given mean value and sigma (distribution's 1st and 2nd order)
2038 % 2) Calculate the expected level of noise from sample statistics
2039 % 3) Calculates the PSD of the individual parts, with:
2040 % Hamming window and linear detrend
2041 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2042 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
2043 % 6) Evaluate the expected value, estimated from the std of the
2044 % individual segments
2045 % 7) Compares with the mean performed on the individual segments
2046 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2047 % 9) Compares the grand-mean with the estimated white noise level
2048
2049 % </TestDescription>
2050 function result = utp_29
2051
2052 % <SyntaxDescription>
2053 %
2054 % 1) Split the reference data into different segments
2055 % 2) Calculates the PSD of the individual parts
2056 %
2057 % </SyntaxDescription>
2058
2059 % <SyntaxCode>
2060 try
2061
2062 % Calculate the expected level of noise from sample statistics
2063
2064 % Evaluate the PSDs on the parts, employing:
2065 % - Hamming window
2066 % - linear detrend
2067 win_type = 'Hamming';
2068 detrend_order = 1;
2069
2070 [S_ap, S_sigma_total, S_expected_total] = ...
2071 test_psd_energy_prepare_data(win_type, detrend_order, plist);
2072
2073 stest = true;
2074
2075 catch err
2076 disp(err.message)
2077 stest = false;
2078 end
2079 % </SyntaxCode>
2080
2081 % <AlgoDescription>
2082 %
2083 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2084 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
2085 % 3) Evaluate the expected value, estimated from the std of the
2086 % individual segments
2087 % 4) Compares with the mean performed on the individual segments
2088 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2089 % 6) Compares the grand-mean with the estimated white noise level
2090
2091 % </AlgoDescription>
2092
2093 % <AlgoCode>
2094
2095 if stest
2096 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
2097 else
2098 atest = false;
2099 end
2100
2101 % </AlgoCode>
2102
2103 % Return a result structure
2104 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2105 end % END UTP_29
2106
2107
2108 %% UTP_30
2109
2110 % <TestDescription>
2111 %
2112 % Tests "conservation of energy":
2113 % 1) white noise produced from normal pdf, with:
2114 % a given mean value and sigma (distribution's 1st and 2nd order)
2115 % 2) Calculate the expected level of noise from sample statistics
2116 % 3) Calculates the PSD of the individual parts, with:
2117 % BH92 window and no detrend
2118 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2119 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
2120 % 6) Evaluate the expected value, estimated from the std of the
2121 % individual segments
2122 % 7) Compares with the mean performed on the individual segments
2123 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2124 % 9) Compares the grand-mean with the estimated white noise level
2125
2126 % </TestDescription>
2127 function result = utp_30
2128
2129 % <SyntaxDescription>
2130 %
2131 % 1) Split the reference data into different segments
2132 % 2) Calculates the PSD of the individual parts
2133 %
2134 % </SyntaxDescription>
2135
2136 % <SyntaxCode>
2137 try
2138
2139 % Calculate the expected level of noise from sample statistics
2140
2141 % Evaluate the PSDs on the parts, employing:
2142 % - BH92 window
2143 % - no detrend
2144 win_type = 'BH92';
2145 detrend_order = -1;
2146
2147 [S_ap, S_sigma_total, S_expected_total] = ...
2148 test_psd_energy_prepare_data(win_type, detrend_order, plist);
2149
2150 stest = true;
2151
2152 catch err
2153 disp(err.message)
2154 stest = false;
2155 end
2156 % </SyntaxCode>
2157
2158 % <AlgoDescription>
2159 %
2160 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2161 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
2162 % 3) Evaluate the expected value, estimated from the std of the
2163 % individual segments
2164 % 4) Compares with the mean performed on the individual segments
2165 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2166 % 6) Compares the grand-mean with the estimated white noise level
2167
2168 % </AlgoDescription>
2169
2170 % <AlgoCode>
2171
2172 if stest
2173 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
2174 else
2175 atest = false;
2176 end
2177
2178 % </AlgoCode>
2179
2180 % Return a result structure
2181 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2182 end % END UTP_30
2183
2184
2185 %% UTP_31
2186
2187 % <TestDescription>
2188 %
2189 % Tests "conservation of energy":
2190 % 1) white noise produced from normal pdf, with:
2191 % a given mean value and sigma (distribution's 1st and 2nd order)
2192 % 2) Calculate the expected level of noise from sample statistics
2193 % 3) Calculates the PSD of the individual parts, with:
2194 % BH92 window and mean detrend
2195 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2196 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
2197 % 6) Evaluate the expected value, estimated from the std of the
2198 % individual segments
2199 % 7) Compares with the mean performed on the individual segments
2200 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2201 % 9) Compares the grand-mean with the estimated white noise level
2202
2203 % </TestDescription>
2204 function result = utp_31
2205
2206 % <SyntaxDescription>
2207 %
2208 % 1) Split the reference data into different segments
2209 % 2) Calculates the PSD of the individual parts
2210 %
2211 % </SyntaxDescription>
2212
2213 % <SyntaxCode>
2214 try
2215
2216 % Calculate the expected level of noise from sample statistics
2217
2218 % Evaluate the PSDs on the parts, employing:
2219 % - BH92 window
2220 % - mean detrend
2221 win_type = 'BH92';
2222 detrend_order = 0;
2223
2224 [S_ap, S_sigma_total, S_expected_total] = ...
2225 test_psd_energy_prepare_data(win_type, detrend_order, plist);
2226
2227 stest = true;
2228
2229 catch err
2230 disp(err.message)
2231 stest = false;
2232 end
2233 % </SyntaxCode>
2234
2235 % <AlgoDescription>
2236 %
2237 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2238 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
2239 % 3) Evaluate the expected value, estimated from the std of the
2240 % individual segments
2241 % 4) Compares with the mean performed on the individual segments
2242 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2243 % 6) Compares the grand-mean with the estimated white noise level
2244
2245 % </AlgoDescription>
2246
2247 % <AlgoCode>
2248
2249 if stest
2250 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
2251 else
2252 atest = false;
2253 end
2254
2255 % </AlgoCode>
2256
2257 % Return a result structure
2258 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2259 end % END UTP_31
2260
2261 %% UTP_32
2262
2263 % <TestDescription>
2264 %
2265 % Tests "conservation of energy":
2266 % 1) white noise produced from normal pdf, with:
2267 % a given mean value and sigma (distribution's 1st and 2nd order)
2268 % 2) Calculate the expected level of noise from sample statistics
2269 % 3) Calculates the PSD of the individual parts, with:
2270 % BH92 window and linear detrend
2271 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2272 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
2273 % 6) Evaluate the expected value, estimated from the std of the
2274 % individual segments
2275 % 7) Compares with the mean performed on the individual segments
2276 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2277 % 9) Compares the grand-mean with the estimated white noise level
2278
2279 % </TestDescription>
2280 function result = utp_32
2281
2282 % <SyntaxDescription>
2283 %
2284 % 1) Split the reference data into different segments
2285 % 2) Calculates the PSD of the individual parts
2286 %
2287 % </SyntaxDescription>
2288
2289 % <SyntaxCode>
2290 try
2291
2292 % Calculate the expected level of noise from sample statistics
2293
2294 % Evaluate the PSDs on the parts, employing:
2295 % - BH92 window
2296 % - linear detrend
2297 win_type = 'BH92';
2298 detrend_order = 1;
2299
2300 [S_ap, S_sigma_total, S_expected_total] = ...
2301 test_psd_energy_prepare_data(win_type, detrend_order, plist);
2302
2303 stest = true;
2304
2305 catch err
2306 disp(err.message)
2307 stest = false;
2308 end
2309 % </SyntaxCode>
2310
2311 % <AlgoDescription>
2312 %
2313 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2314 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
2315 % 3) Evaluate the expected value, estimated from the std of the
2316 % individual segments
2317 % 4) Compares with the mean performed on the individual segments
2318 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2319 % 6) Compares the grand-mean with the estimated white noise level
2320
2321 % </AlgoDescription>
2322
2323 % <AlgoCode>
2324
2325 if stest
2326 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
2327 else
2328 atest = false;
2329 end
2330
2331 % </AlgoCode>
2332
2333 % Return a result structure
2334 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2335 end % END UTP_32
2336
2337 %% UTP_33
2338
2339 % <TestDescription>
2340 %
2341 % Tests "conservation of energy":
2342 % 1) white noise produced from normal pdf, with:
2343 % a given mean value and sigma (distribution's 1st and 2nd order)
2344 % 2) Calculate the expected level of noise from sample statistics
2345 % 3) Calculates the PSD of the individual parts, with:
2346 % Kaiser200 window and no detrend
2347 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2348 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
2349 % 6) Evaluate the expected value, estimated from the std of the
2350 % individual segments
2351 % 7) Compares with the mean performed on the individual segments
2352 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2353 % 9) Compares the grand-mean with the estimated white noise level
2354
2355 % </TestDescription>
2356 function result = utp_33
2357
2358 % <SyntaxDescription>
2359 %
2360 % 1) Split the reference data into different segments
2361 % 2) Calculates the PSD of the individual parts
2362 %
2363 % </SyntaxDescription>
2364
2365 % <SyntaxCode>
2366 try
2367
2368 % Calculate the expected level of noise from sample statistics
2369
2370 % Evaluate the PSDs on the parts, employing:
2371 % - Kaiser200 window
2372 % - no detrend
2373 win_type = 'Kaiser';
2374 psll = 200;
2375 detrend_order = -1;
2376
2377 [S_ap, S_sigma_total, S_expected_total] = ...
2378 test_psd_energy_prepare_data(win_type, detrend_order, plist('psll', psll));
2379
2380 stest = true;
2381
2382 catch err
2383 disp(err.message)
2384 stest = false;
2385 end
2386 % </SyntaxCode>
2387
2388 % <AlgoDescription>
2389 %
2390 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2391 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
2392 % 3) Evaluate the expected value, estimated from the std of the
2393 % individual segments
2394 % 4) Compares with the mean performed on the individual segments
2395 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2396 % 6) Compares the grand-mean with the estimated white noise level
2397
2398 % </AlgoDescription>
2399
2400 % <AlgoCode>
2401
2402 if stest
2403 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
2404 else
2405 atest = false;
2406 end
2407
2408 % </AlgoCode>
2409
2410 % Return a result structure
2411 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2412 end % END UTP_33
2413
2414
2415 %% UTP_34
2416
2417 % <TestDescription>
2418 %
2419 % Tests "conservation of energy":
2420 % 1) white noise produced from normal pdf, with:
2421 % a given mean value and sigma (distribution's 1st and 2nd order)
2422 % 2) Calculate the expected level of noise from sample statistics
2423 % 3) Calculates the PSD of the individual parts, with:
2424 % Kaiser200 window and mean detrend
2425 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2426 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
2427 % 6) Evaluate the expected value, estimated from the std of the
2428 % individual segments
2429 % 7) Compares with the mean performed on the individual segments
2430 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2431 % 9) Compares the grand-mean with the estimated white noise level
2432
2433 % </TestDescription>
2434 function result = utp_34
2435
2436 % <SyntaxDescription>
2437 %
2438 % 1) Split the reference data into different segments
2439 % 2) Calculates the PSD of the individual parts
2440 %
2441 % </SyntaxDescription>
2442
2443 % <SyntaxCode>
2444 try
2445
2446 % Calculate the expected level of noise from sample statistics
2447
2448 % Evaluate the PSDs on the parts, employing:
2449 % - Kaiser200 window
2450 % - mean detrend
2451 win_type = 'Kaiser';
2452 psll = 200;
2453 detrend_order = 0;
2454
2455 [S_ap, S_sigma_total, S_expected_total] = ...
2456 test_psd_energy_prepare_data(win_type, detrend_order, plist('psll', psll));
2457
2458 stest = true;
2459
2460 catch err
2461 disp(err.message)
2462 stest = false;
2463 end
2464 % </SyntaxCode>
2465
2466 % <AlgoDescription>
2467 %
2468 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2469 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
2470 % 3) Evaluate the expected value, estimated from the std of the
2471 % individual segments
2472 % 4) Compares with the mean performed on the individual segments
2473 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2474 % 6) Compares the grand-mean with the estimated white noise level
2475
2476 % </AlgoDescription>
2477
2478 % <AlgoCode>
2479
2480 if stest
2481 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
2482 else
2483 atest = false;
2484 end
2485
2486 % </AlgoCode>
2487
2488 % Return a result structure
2489 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2490 end % END UTP_34
2491
2492 %% UTP_35
2493
2494 % <TestDescription>
2495 %
2496 % Tests "conservation of energy":
2497 % 1) white noise produced from normal pdf, with:
2498 % a given mean value and sigma (distribution's 1st and 2nd order)
2499 % 2) Calculate the expected level of noise from sample statistics
2500 % 3) Calculates the PSD of the individual parts, with:
2501 % Kaiser200 window and linear detrend
2502 % 4) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2503 % 5) Checks on individual PSDs: mean and standard deviation of the PSD points
2504 % 6) Evaluate the expected value, estimated from the std of the
2505 % individual segments
2506 % 7) Compares with the mean performed on the individual segments
2507 % 8) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2508 % 9) Compares the grand-mean with the estimated white noise level
2509
2510 % </TestDescription>
2511 function result = utp_35
2512
2513 % <SyntaxDescription>
2514 %
2515 % 1) Split the reference data into different segments
2516 % 2) Calculates the PSD of the individual parts
2517 %
2518 % </SyntaxDescription>
2519
2520 % <SyntaxCode>
2521 try
2522
2523 % Calculate the expected level of noise from sample statistics
2524
2525 % Evaluate the PSDs on the parts, employing:
2526 % - Kaiser200 window
2527 % - linear detrend
2528 win_type = 'Kaiser';
2529 psll = 200;
2530 detrend_order = 1;
2531
2532 [S_ap, S_sigma_total, S_expected_total] = ...
2533 test_psd_energy_prepare_data(win_type, detrend_order, plist('psll', psll));
2534
2535 stest = true;
2536
2537 catch err
2538 disp(err.message)
2539 stest = false;
2540 end
2541 % </SyntaxCode>
2542
2543 % <AlgoDescription>
2544 %
2545 % 1) Evaluate the mean and standard deviation of the Checks on individual PSDs:
2546 % 2) Checks on individual PSDs: mean and standard deviation of the PSD points
2547 % 3) Evaluate the expected value, estimated from the std of the
2548 % individual segments
2549 % 4) Compares with the mean performed on the individual segments
2550 % 5) Checks on averaged PSDs: mean and standard deviation of all the PSD points
2551 % 6) Compares the grand-mean with the estimated white noise level
2552
2553 % </AlgoDescription>
2554
2555 % <AlgoCode>
2556
2557 if stest
2558 atest = algo_test_psd_energy(ap, S_ap, S_expected_total, plist);
2559 else
2560 atest = false;
2561 end
2562
2563 % </AlgoCode>
2564
2565 % Return a result structure
2566 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2567 end % END UTP_35
2568
2569 %% UTP_38
2570
2571 % <TestDescription>
2572 %
2573 % Tests "conservation of energy":
2574 % 1) white noise produced from normal pdf, with:
2575 % a given mean value and sigma (distribution's 1st and 2nd order)
2576 % 2) evaluate the sample mean value m and standard deviation s
2577 % 3) add a given trend of order n
2578 % 4) psd of the noise, with proper detrending
2579 % 5) compares the sqrt(sum(S * df)) with the standard deviation s
2580 %
2581
2582 % </TestDescription>
2583 function result = utp_38
2584
2585 % <SyntaxDescription>
2586 %
2587 % 1) Prepare the test tsdata:
2588 % white noise from normal distribution + offset
2589 % 2) Calculate the statistical parameters
2590 % 3) Add a trend
2591 % 4) Estimate the psd
2592 %
2593 % </SyntaxDescription>
2594
2595 % <SyntaxCode>
2596 try
2597 fs = 10;
2598 nsecs = 3600;
2599 sigma_distr = 1;
2600 trend_0 = 0;
2601 trend_1 = 2e-6;
2602 trend_2 = 5e-10;
2603
2604 % White noise
2605 type = 'Normal';
2606 a_n = ao(plist('waveform', 'noise', ...
2607 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
2608
2609 % Drift signal
2610 a_d = ao(plist('tsfcn', [num2str(trend_0) '*t.^0 + ' num2str(trend_1) ' *t + ' num2str(trend_2) ' *t.^2'], ...
2611 'fs', fs, 'nsecs', nsecs));
2612
2613 % Total signal
2614 a = a_n + a_d;
2615
2616 % Estimate the mean value and the sigma of the white noise time-series data
2617 sigma_sample = std(a_n.y, 1);
2618
2619 % Evaluate the psd of the white noise time-series data
2620 % Parameters are crucial for the estimate to be correct!
2621 win = 'Rectangular';
2622 olap = 0;
2623 detrend = 0;
2624 n_pts = -1;
2625 scale = 'PSD';
2626
2627 S_n = a_n.psd(plist('Win', win, 'olap', olap, ...
2628 'Nfft', n_pts, 'order', detrend, 'scale', scale));
2629
2630 % Evaluate the psd of the "drifting" time-series data
2631 % Parameters are crucial for the estimate to be correct!
2632 win = 'Rectangular';
2633 olap = 0;
2634 detrend = 2;
2635 n_pts = -1;
2636 scale = 'PSD';
2637
2638 S = a.psd(plist('Win', win, 'olap', olap, ...
2639 'Nfft', n_pts, 'order', detrend, 'scale', scale));
2640 stest = true;
2641
2642 catch err
2643 disp(err.message)
2644 stest = false;
2645 end
2646 % </SyntaxCode>
2647
2648 % <AlgoDescription>
2649 %
2650 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
2651 %
2652 % </AlgoDescription>
2653
2654 % <AlgoCode>
2655 atest = true;
2656 TOL = 1e-3;
2657
2658 if stest
2659 % Integrals of the spectra
2660 sigma_psd = sqrt(sum(S.y * S.x(2)));
2661 sigma_psd_n = sqrt(sum(S_n.y * S_n.x(2)));
2662
2663 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL || ...
2664 abs(sigma_psd_n - sigma_sample) / mean([sigma_psd_n sigma_sample]) >= TOL
2665 atest = false;
2666 end
2667 else
2668 atest = false;
2669 end
2670 % </AlgoCode>
2671
2672 % Return a result structure
2673 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2674 end % END UTP_38
2675
2676
2677 %% UTP_39
2678
2679 % <TestDescription>
2680 %
2681 % Tests "conservation of energy":
2682 % 1) white noise produced from uniform pdf, with:
2683 % a given mean value and sigma (distribution's 1st and 2nd order)
2684 % 2) evaluate the sample mean value m and standard deviation s
2685 % 3) add a given trend of order n
2686 % 4) psd of the noise, with proper detrending
2687 % 5) compares the sqrt(sum(S * df)) with the standard deviation s
2688 %
2689
2690 % </TestDescription>
2691 function result = utp_39
2692
2693 % <SyntaxDescription>
2694 %
2695 % 1) Prepare the test tsdata:
2696 % white noise from uniform distribution + offset
2697 % 2) Calculate the statistical parameters
2698 % 3) Add a trend
2699 % 4) Estimate the psd
2700 %
2701 % </SyntaxDescription>
2702
2703 % <SyntaxCode>
2704 try
2705 % Build time-series test data
2706 fs = 1;
2707 nsecs = 86400;
2708 sigma_distr = 1e-9;
2709 trend_0 = 3e-7;
2710 trend_1 = 1e-7;
2711 trend_2 = 2e-14;
2712
2713 % White noise
2714 type = 'Uniform';
2715 a_n = ao(plist('waveform', 'noise', ...
2716 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
2717
2718 % Drift signal
2719 a_d = ao(plist('tsfcn', [num2str(trend_0) '*t.^0 + ' num2str(trend_1) ' *t + ' num2str(trend_2) ' *t.^2'], ...
2720 'fs', fs, 'nsecs', nsecs));
2721
2722 % Total signal
2723 a = a_n + a_d;
2724
2725 % Estimate the mean value and the sigma of the white noise time-series data
2726 sigma_sample = std(a_n.y, 1);
2727
2728 % Evaluate the psd of the white noise time-series data
2729 % Parameters are crucial for the estimate to be correct!
2730 win = 'Rectangular';
2731 olap = 0;
2732 detrend = 0;
2733 n_pts = -1;
2734 scale = 'PSD';
2735
2736 S_n = a_n.psd(plist('Win', win, 'olap', olap, ...
2737 'Nfft', n_pts, 'order', detrend, 'scale', scale));
2738
2739 % Evaluate the psd of the "drifting" time-series data
2740 % Parameters are crucial for the estimate to be correct!
2741 win = 'Rectangular';
2742 olap = 0;
2743 detrend = 2;
2744 n_pts = -1;
2745 scale = 'PSD';
2746
2747 S = a.psd(plist('Win', win, 'olap', olap, ...
2748 'Nfft', n_pts, 'order', detrend, 'scale', scale));
2749 stest = true;
2750
2751 catch err
2752 disp(err.message)
2753 stest = false;
2754 end
2755 % </SyntaxCode>
2756
2757 % <AlgoDescription>
2758 %
2759 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
2760 %
2761 % </AlgoDescription>
2762
2763 % <AlgoCode>
2764 atest = true;
2765 TOL = 1e-3;
2766
2767 if stest
2768 % Integrals of the spectra
2769 sigma_psd = sqrt(sum(S.y * S.x(2)));
2770 sigma_psd_n = sqrt(sum(S_n.y * S_n.x(2)));
2771
2772 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL || ...
2773 abs(sigma_psd_n - sigma_sample) / mean([sigma_psd_n sigma_sample]) >= TOL
2774 atest = false;
2775 end
2776 else
2777 atest = false;
2778 end
2779 % </AlgoCode>
2780
2781 % Return a result structure
2782 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2783 end % END UTP_39
2784
2785 %% UTP_40
2786
2787 % <TestDescription>
2788 %
2789 % Tests "conservation of energy":
2790 % 1) white noise produced from normal pdf, with:
2791 % a given mean value and sigma (distribution's 1st and 2nd order)
2792 % 2) evaluate the sample mean value m and standard deviation s
2793 % 3) add a given trend of order n
2794 % 4) psd of the noise, with proper detrending
2795 % 5) compares the sqrt(sum(S * df)) with the standard deviation s
2796 %
2797
2798 % </TestDescription>
2799 function result = utp_40
2800
2801 % <SyntaxDescription>
2802 %
2803 % 1) Prepare the test tsdata:
2804 % white noise from normal distribution + offset
2805 % 2) Calculate the statistical parameters
2806 % 3) Add a trend
2807 % 4) Estimate the psd
2808 %
2809 % </SyntaxDescription>
2810
2811 % <SyntaxCode>
2812 try
2813 % Build time-series test data
2814
2815 fs = 100;
2816 nsecs = 100;
2817 sigma_distr = 1e-12;
2818 trend_0 = -5e-12;
2819 trend_1 = 4e-13;
2820 trend_2 = 1.17e-14;
2821
2822 % White noise
2823 type = 'Normal';
2824 a_n = ao(plist('waveform', 'noise', ...
2825 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
2826
2827 % Drift signal
2828 a_d = ao(plist('tsfcn', [num2str(trend_0) '*t.^0 + ' num2str(trend_1) ' *t + ' num2str(trend_2) ' *t.^2'], ...
2829 'fs', fs, 'nsecs', nsecs));
2830
2831 % Total signal
2832 a = a_n + a_d;
2833
2834 % Estimate the mean value and the sigma of the white noise time-series data
2835 sigma_sample = std(a_n.y, 1);
2836
2837 % Evaluate the psd of the white noise time-series data
2838 % Parameters are crucial for the estimate to be correct!
2839 win = 'Rectangular';
2840 olap = 0;
2841 detrend = 0;
2842 n_pts = -1;
2843 scale = 'PSD';
2844
2845 S_n = a_n.psd(plist('Win', win, 'olap', olap, ...
2846 'Nfft', n_pts, 'order', detrend, 'scale', scale));
2847
2848 % Evaluate the psd of the "drifting" time-series data
2849 % Parameters are crucial for the estimate to be correct!
2850 win = 'Rectangular';
2851 olap = 0;
2852 detrend = 2;
2853 n_pts = -1;
2854 scale = 'PSD';
2855
2856 S = a.psd(plist('Win', win, 'olap', olap, ...
2857 'Nfft', n_pts, 'order', detrend, 'scale', scale));
2858 stest = true;
2859
2860 catch err
2861 disp(err.message)
2862 stest = false;
2863 end
2864 % </SyntaxCode>
2865
2866 % <AlgoDescription>
2867 %
2868 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
2869 %
2870 % </AlgoDescription>
2871
2872 % <AlgoCode>
2873 atest = true;
2874 TOL = 1e-3;
2875
2876 if stest
2877 % Integrals of the spectra
2878 sigma_psd = sqrt(sum(S.y * S.x(2)));
2879 sigma_psd_n = sqrt(sum(S_n.y * S_n.x(2)));
2880
2881 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL || ...
2882 abs(sigma_psd_n - sigma_sample) / mean([sigma_psd_n sigma_sample]) >= TOL
2883 atest = false;
2884 end
2885 else
2886 atest = false;
2887 end
2888 % </AlgoCode>
2889
2890 % Return a result structure
2891 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2892 end % END UTP_40
2893
2894
2895
2896 %% UTP_41
2897
2898 % <TestDescription>
2899 %
2900 % Tests the effect of windowing:
2901 % 1) white noise produced from normal pdf, with:
2902 % a given mean value and sigma (distribution's 1st and 2nd order)
2903 % 2) psd of the noise, without detrending, Rectangular window
2904 % 3) Apply the detrending and the window manually
2905 % 4) psd of the noise, without detrending, Rectangular window
2906 % 5) compares the to psds
2907 %
2908
2909 % </TestDescription>
2910 function result = utp_41
2911
2912 % <SyntaxDescription>
2913 %
2914 % 1) Prepare the test tsdata:
2915 % white noise from normal distribution + offset
2916 % 2) Calculate the statistical parameters
2917 % 3) Estimate the psd without detrending, Rectangular window
2918 % 4) Manually apply window to the data
2919 % 5) Estimate the psd without detrending, Rectangular window
2920 %
2921 % </SyntaxDescription>
2922
2923 % <SyntaxCode>
2924 try
2925
2926 noise_type = 'Normal';
2927 win_type = 'Rectangular';
2928 scale = 'PSD';
2929
2930 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
2931 plist('win_test', true));
2932
2933 stest = true;
2934
2935 catch err
2936 disp(err.message)
2937 stest = false;
2938 end
2939 % </SyntaxCode>
2940
2941 % <AlgoDescription>
2942 %
2943 % 1) Check that calculated psds are identical within a part in 10^12
2944 %
2945 % </AlgoDescription>
2946
2947 % <AlgoCode>
2948
2949 TOL = 1e-12;
2950
2951 if stest
2952 atest = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
2953 else
2954 atest = false;
2955 end
2956
2957 % </AlgoCode>
2958
2959 % Return a result structure
2960 result = utp_prepare_result(atest, stest, dbstack, mfilename);
2961 end % END UTP_41
2962
2963 %% UTP_42
2964
2965 % <TestDescription>
2966 %
2967 % Tests the effect of windowing:
2968 % 1) white noise produced from normal pdf, with:
2969 % a given mean value and sigma (distribution's 1st and 2nd order)
2970 % 2) psd of the noise, without detrending, BH92 window
2971 % 3) Apply the detrending and the window manually
2972 % 4) psd of the noise, without detrending, Rectangular window
2973 % 5) compares the to psds
2974 %
2975
2976 % </TestDescription>
2977 function result = utp_42
2978
2979 % <SyntaxDescription>
2980 %
2981 % 1) Prepare the test tsdata:
2982 % white noise from normal distribution + offset
2983 % 2) Calculate the statistical parameters
2984 % 3) Estimate the psd without detrending, BH92 window
2985 % 4) Manually apply window to the data
2986 % 5) Estimate the psd without detrending, Rectangular window
2987 %
2988 % </SyntaxDescription>
2989
2990 % <SyntaxCode>
2991 try
2992
2993 noise_type = 'Normal';
2994 win_type = 'BH92';
2995 scale = 'PSD';
2996
2997 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
2998 plist('win_test', true));
2999
3000 stest = true;
3001
3002 catch err
3003 disp(err.message)
3004 stest = false;
3005 end
3006 % </SyntaxCode>
3007
3008 % <AlgoDescription>
3009 %
3010 % 1) Check that calculated psds are identical within a part in 10^12
3011 %
3012 % </AlgoDescription>
3013
3014 % <AlgoCode>
3015
3016 TOL = 1e-12;
3017
3018 if stest
3019 atest = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3020 else
3021 atest = false;
3022 end
3023
3024 % </AlgoCode>
3025
3026 % Return a result structure
3027 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3028 end % END UTP_42
3029
3030
3031 %% UTP_43
3032
3033 % <TestDescription>
3034 %
3035 % Tests the effect of windowing:
3036 % 1) white noise produced from normal pdf, with:
3037 % a given mean value and sigma (distribution's 1st and 2nd order)
3038 % 2) psd of the noise, without detrending, Hamming window
3039 % 3) Apply the detrending and the window manually
3040 % 4) psd of the noise, without detrending, Rectangular window
3041 % 5) compares the to psds
3042 %
3043
3044 % </TestDescription>
3045 function result = utp_43
3046
3047 % <SyntaxDescription>
3048 %
3049 % 1) Prepare the test tsdata:
3050 % white noise from normal distribution + offset
3051 % 2) Calculate the statistical parameters
3052 % 3) Estimate the psd without detrending, Hamming window
3053 % 4) Manually apply window to the data
3054 % 5) Estimate the psd without detrending, Rectangular window
3055 %
3056 % </SyntaxDescription>
3057
3058 % <SyntaxCode>
3059 try
3060
3061 noise_type = 'Normal';
3062 win_type = 'Hamming';
3063 scale = 'PSD';
3064
3065 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3066 plist('win_test', true));
3067
3068 stest = true;
3069
3070 catch err
3071 disp(err.message)
3072 stest = false;
3073 end
3074 % </SyntaxCode>
3075
3076 % <AlgoDescription>
3077 %
3078 % 1) Check that calculated psds are identical within a part in 10^12
3079 %
3080 % </AlgoDescription>
3081
3082 % <AlgoCode>
3083
3084 TOL = 1e-12;
3085
3086 if stest
3087 atest = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3088 else
3089 atest = false;
3090 end
3091
3092 % </AlgoCode>
3093
3094 % Return a result structure
3095 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3096 end % END UTP_43
3097
3098 %% UTP_44
3099
3100 % <TestDescription>
3101 %
3102 % Tests the effect of windowing:
3103 % 1) white noise produced from normal pdf, with:
3104 % a given mean value and sigma (distribution's 1st and 2nd order)
3105 % 2) psd of the noise, without detrending, Hanning window
3106 % 3) Apply the detrending and the window manually
3107 % 4) psd of the noise, without detrending, Rectangular window
3108 % 5) compares the to psds
3109 %
3110
3111 % </TestDescription>
3112 function result = utp_44
3113
3114 % <SyntaxDescription>
3115 %
3116 % 1) Prepare the test tsdata:
3117 % white noise from normal distribution + offset
3118 % 2) Calculate the statistical parameters
3119 % 3) Estimate the psd without detrending, Hanning window
3120 % 4) Manually apply window to the data
3121 % 5) Estimate the psd without detrending, Rectangular window
3122 %
3123 % </SyntaxDescription>
3124
3125 % <SyntaxCode>
3126 msg = '';
3127 try
3128
3129 noise_type = 'Normal';
3130 win_type = 'Hanning';
3131 scale = 'PSD';
3132
3133 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3134 plist('win_test', true));
3135
3136 stest = true;
3137
3138 catch err
3139 disp(err.message)
3140 msg = err.message;
3141 stest = false;
3142 end
3143 % </SyntaxCode>
3144
3145 % <AlgoDescription>
3146 %
3147 % 1) Check that calculated psds are identical within a part in 10^12
3148 %
3149 % </AlgoDescription>
3150
3151 % <AlgoCode>
3152
3153 TOL = 1e-12;
3154
3155 if stest
3156 [atest, msg] = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3157 else
3158 atest = false;
3159 end
3160
3161 % </AlgoCode>
3162
3163 % Return a result structure
3164 result = utp_prepare_result(atest, stest, dbstack, mfilename, msg);
3165 end % END UTP_44
3166
3167 %% UTP_45
3168
3169 % <TestDescription>
3170 %
3171 % Tests the effect of windowing:
3172 % 1) white noise produced from normal pdf, with:
3173 % a given mean value and sigma (distribution's 1st and 2nd order)
3174 % 2) psd of the noise, without detrending, Bartlett window
3175 % 3) Apply the detrending and the window manually
3176 % 4) psd of the noise, without detrending, Rectangular window
3177 % 5) compares the to psds
3178 %
3179
3180 % </TestDescription>
3181 function result = utp_45
3182
3183 % <SyntaxDescription>
3184 %
3185 % 1) Prepare the test tsdata:
3186 % white noise from normal distribution + offset
3187 % 2) Calculate the statistical parameters
3188 % 3) Estimate the psd without detrending, Bartlett window
3189 % 4) Manually apply window to the data
3190 % 5) Estimate the psd without detrending, Rectangular window
3191 %
3192 % </SyntaxDescription>
3193
3194 % <SyntaxCode>
3195 try
3196
3197 noise_type = 'Normal';
3198 win_type = 'Bartlett';
3199 scale = 'PSD';
3200
3201 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3202 plist('win_test', true));
3203
3204 stest = true;
3205
3206 catch err
3207 disp(err.message)
3208 stest = false;
3209 end
3210 % </SyntaxCode>
3211
3212 % <AlgoDescription>
3213 %
3214 % 1) Check that calculated psds are identical within a part in 10^12
3215 %
3216 % </AlgoDescription>
3217
3218 % <AlgoCode>
3219
3220 TOL = 1e-12;
3221
3222 if stest
3223 atest = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3224 else
3225 atest = false;
3226 end
3227
3228 % </AlgoCode>
3229
3230 % Return a result structure
3231 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3232 end % END UTP_45
3233
3234 %% UTP_46
3235
3236 % <TestDescription>
3237 %
3238 % Tests the effect of windowing:
3239 % 1) white noise produced from normal pdf, with:
3240 % a given mean value and sigma (distribution's 1st and 2nd order)
3241 % 2) psd of the noise, without detrending, Nuttall3 window
3242 % 3) Apply the detrending and the window manually
3243 % 4) psd of the noise, without detrending, Rectangular window
3244 % 5) compares the to psds
3245 %
3246
3247 % </TestDescription>
3248 function result = utp_46
3249
3250 % <SyntaxDescription>
3251 %
3252 % 1) Prepare the test tsdata:
3253 % white noise from normal distribution + offset
3254 % 2) Calculate the statistical parameters
3255 % 3) Estimate the psd without detrending, Nuttall3 window
3256 % 4) Manually apply window to the data
3257 % 5) Estimate the psd without detrending, Rectangular window
3258 %
3259 % </SyntaxDescription>
3260
3261 % <SyntaxCode>
3262 msg = '';
3263 try
3264
3265 noise_type = 'Normal';
3266 win_type = 'Nuttall3';
3267 scale = 'PSD';
3268
3269 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3270 plist('win_test', true));
3271
3272 stest = true;
3273
3274 catch err
3275 disp(err.message)
3276 msg = err.message;
3277 stest = false;
3278 end
3279 % </SyntaxCode>
3280
3281 % <AlgoDescription>
3282 %
3283 % 1) Check that calculated psds are identical within a part in 10^12
3284 %
3285 % </AlgoDescription>
3286
3287 % <AlgoCode>
3288
3289 TOL = 1e-12;
3290
3291 if stest
3292 [atest, msg] = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3293 else
3294 atest = false;
3295 end
3296
3297 % </AlgoCode>
3298
3299 % Return a result structure
3300 result = utp_prepare_result(atest, stest, dbstack, mfilename, msg);
3301 end % END UTP_46
3302
3303 %% UTP_47
3304
3305 % <TestDescription>
3306 %
3307 % Tests the effect of windowing:
3308 % 1) white noise produced from normal pdf, with:
3309 % a given mean value and sigma (distribution's 1st and 2nd order)
3310 % 2) psd of the noise, without detrending, Kaiser psll = random window
3311 % 3) Apply the detrending and the window manually
3312 % 4) psd of the noise, without detrending, Rectangular window
3313 % 5) compares the to psds
3314 %
3315
3316 % </TestDescription>
3317 function result = utp_47
3318
3319 % <SyntaxDescription>
3320 %
3321 % 1) Prepare the test tsdata:
3322 % white noise from normal distribution + offset
3323 % 2) Calculate the statistical parameters
3324 % 3) Estimate the psd without detrending, Kaiser psll = random window
3325 % 4) Manually apply window to the data
3326 % 5) Estimate the psd without detrending, Rectangular window
3327 %
3328 % </SyntaxDescription>
3329
3330 % <SyntaxCode>
3331 try
3332
3333 noise_type = 'Normal';
3334 win_type = 'Kaiser';
3335 psll = utils.math.randelement([10:10:200],1);
3336 scale = 'PSD';
3337
3338 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3339 plist('win_test', true));
3340
3341 stest = true;
3342
3343 catch err
3344 disp(err.message)
3345 stest = false;
3346 end
3347 % </SyntaxCode>
3348
3349 % <AlgoDescription>
3350 %
3351 % 1) Check that calculated psds are identical within a part in 10^12
3352 %
3353 % </AlgoDescription>
3354
3355 % <AlgoCode>
3356
3357 TOL = 1e-12;
3358
3359 if stest
3360 atest = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3361 else
3362 atest = false;
3363 end
3364
3365 % </AlgoCode>
3366
3367 % Return a result structure
3368 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3369 end % END UTP_47
3370
3371 %% UTP_48
3372
3373 % <TestDescription>
3374 %
3375 % Tests the effect of windowing:
3376 % 1) white noise produced from normal pdf, with:
3377 % a given mean value and sigma (distribution's 1st and 2nd order)
3378 % 2) psd of the noise, without detrending, Kaiser psll = default window
3379 % 3) Apply the detrending and the window manually
3380 % 4) psd of the noise, without detrending, Rectangular window
3381 % 5) compares the to psds
3382 %
3383
3384 % </TestDescription>
3385 function result = utp_48
3386
3387 % <SyntaxDescription>
3388 %
3389 % 1) Prepare the test tsdata:
3390 % white noise from normal distribution + offset
3391 % 2) Calculate the statistical parameters
3392 % 3) Estimate the psd without detrending, Kaiser psll = default window
3393 % 4) Manually apply window to the data
3394 % 5) Estimate the psd without detrending, Rectangular window
3395 %
3396 % </SyntaxDescription>
3397
3398 % <SyntaxCode>
3399 try
3400
3401 noise_type = 'Normal';
3402 win_type = 'Kaiser';
3403 scale = 'PSD';
3404
3405 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3406 plist('win_test', true));
3407
3408 stest = true;
3409
3410 catch err
3411 disp(err.message)
3412 stest = false;
3413 end
3414 % </SyntaxCode>
3415
3416 % <AlgoDescription>
3417 %
3418 % 1) Check that calculated psds are identical within a part in 10^12
3419 %
3420 % </AlgoDescription>
3421
3422 % <AlgoCode>
3423
3424 TOL = 1e-12;
3425
3426 if stest
3427 atest = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3428 else
3429 atest = false;
3430 end
3431
3432 % </AlgoCode>
3433
3434 % Return a result structure
3435 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3436 end % END UTP_48
3437
3438
3439 %% UTP_49
3440
3441 % <TestDescription>
3442 %
3443 % Tests the effect of windowing:
3444 % 1) white noise produced from normal pdf, with:
3445 % a given mean value and sigma (distribution's 1st and 2nd order)
3446 % 2) psd of the noise, without detrending, Nuttall4 window
3447 % 3) Apply the detrending and the window manually
3448 % 4) psd of the noise, without detrending, Rectangular window
3449 % 5) compares the to psds
3450 %
3451
3452 % </TestDescription>
3453 function result = utp_49
3454
3455 % <SyntaxDescription>
3456 %
3457 % 1) Prepare the test tsdata:
3458 % white noise from normal distribution + offset
3459 % 2) Calculate the statistical parameters
3460 % 3) Estimate the psd without detrending, Nuttall4 window
3461 % 4) Manually apply window to the data
3462 % 5) Estimate the psd without detrending, Rectangular window
3463 %
3464 % </SyntaxDescription>
3465
3466 % <SyntaxCode>
3467 try
3468
3469 noise_type = 'Normal';
3470 win_type = 'Nuttall4';
3471 scale = 'PSD';
3472
3473 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3474 plist('win_test', true));
3475
3476 stest = true;
3477
3478 catch err
3479 disp(err.message)
3480 stest = false;
3481 end
3482 % </SyntaxCode>
3483
3484 % <AlgoDescription>
3485 %
3486 % 1) Check that calculated psds are identical within a part in 10^12
3487 %
3488 % </AlgoDescription>
3489
3490 % <AlgoCode>
3491
3492 TOL = 1e-12;
3493
3494 if stest
3495 atest = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3496 else
3497 atest = false;
3498 end
3499
3500 % </AlgoCode>
3501
3502 % Return a result structure
3503 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3504 end % END UTP_49
3505
3506 %% UTP_50
3507
3508 % <TestDescription>
3509 %
3510 % Tests the effect of windowing:
3511 % 1) white noise produced from normal pdf, with:
3512 % a given mean value and sigma (distribution's 1st and 2nd order)
3513 % 2) psd of the noise, without detrending, SFT3F window
3514 % 3) Apply the detrending and the window manually
3515 % 4) psd of the noise, without detrending, Rectangular window
3516 % 5) compares the to psds
3517 %
3518
3519 % </TestDescription>
3520 function result = utp_50
3521
3522 % <SyntaxDescription>
3523 %
3524 % 1) Prepare the test tsdata:
3525 % white noise from normal distribution + offset
3526 % 2) Calculate the statistical parameters
3527 % 3) Estimate the psd without detrending, SFT3F window
3528 % 4) Manually apply window to the data
3529 % 5) Estimate the psd without detrending, Rectangular window
3530 %
3531 % </SyntaxDescription>
3532
3533 % <SyntaxCode>
3534 try
3535
3536 noise_type = 'Normal';
3537 win_type = 'SFT3F';
3538 scale = 'PSD';
3539
3540 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3541 plist('win_test', true));
3542
3543 stest = true;
3544
3545 catch err
3546 disp(err.message)
3547 stest = false;
3548 end
3549 % </SyntaxCode>
3550
3551 % <AlgoDescription>
3552 %
3553 % 1) Check that calculated psds are identical within a part in 10^12
3554 %
3555 % </AlgoDescription>
3556
3557 % <AlgoCode>
3558
3559 TOL = 1e-12;
3560
3561 if stest
3562 atest = algo_test_psd_win(S_name, S_obj, S_man, plist('TOL', TOL));
3563 else
3564 atest = false;
3565 end
3566
3567 % </AlgoCode>
3568
3569 % Return a result structure
3570 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3571 end % END UTP_50
3572
3573 %% UTP_51
3574
3575 % <TestDescription>
3576 %
3577 % Tests the possibility to set the number of averages rather than setting the Nfft:
3578 % 1) white noise produced from normal pdf, with:
3579 % a given mean value and sigma (distribution's 1st and 2nd order)
3580 % 2) psd of the noise, without detrending, random window, set number of
3581 % averages
3582 % 3) check the effective number of averages
3583 %
3584
3585 % </TestDescription>
3586 function result = utp_51
3587
3588 % <SyntaxDescription>
3589 %
3590 % 1) Prepare the test tsdata:
3591 % white noise from normal distribution + offset
3592 % 2) psd of the noise, without detrending, random window, set number of
3593 % averages
3594 %
3595 % </SyntaxDescription>
3596
3597 % <SyntaxCode>
3598 try
3599 noise_type = 'Normal';
3600
3601 % Evaluate the psd of the white noise time-series data
3602 [a, S1, S2, navs] = prepare_analyze_noise_navs(noise_type, plist);
3603
3604 stest = true;
3605
3606 catch err
3607 disp(err.message)
3608 stest = false;
3609 end
3610 % </SyntaxCode>
3611
3612 % <AlgoDescription>
3613 %
3614 % 1) Check that calculated navs are identical to those requested
3615 %
3616 % </AlgoDescription>
3617
3618 % <AlgoCode>
3619 atest = true;
3620
3621 if stest
3622 % Compare the navs written in the output object with the requested one
3623 if ne(navs, S1.data.navs)
3624 if ne(find(S1.hist.plistUsed, 'navs'), S1.data.navs)
3625 atest = false;
3626 end
3627 end
3628 else
3629 atest = false;
3630 end
3631 % </AlgoCode>
3632
3633 % Return a result structure
3634 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3635 end % END UTP_51
3636
3637 %% UTP_52
3638
3639 % <TestDescription>
3640 %
3641 % Tests the possibility to set the number of averages rather than setting the Nfft:
3642 % 1) white noise produced from normal pdf, with:
3643 % a given mean value and sigma (distribution's 1st and 2nd order)
3644 % 2) psd of the noise, without detrending, random window, random navs
3645 % 3) get the number of averages
3646 % 4) get the nfft used
3647 % 5) run psd again, with the nfft used
3648 % 6) compare the calculated objects
3649 %
3650
3651 % </TestDescription>
3652 function result = utp_52
3653
3654 % <SyntaxDescription>
3655 %
3656 % 1) white noise produced from normal pdf, with:
3657 % a given mean value and sigma (distribution's 1st and 2nd order)
3658 % 2) psd of the noise, without detrending, random window, random navs
3659 % 3) get the number of averages
3660 % 4) get the nfft used
3661 % 5) run psd again, with the nfft used
3662 %
3663 % </SyntaxDescription>
3664
3665 % <SyntaxCode>
3666 try
3667 noise_type = 'Normal';
3668
3669 % Evaluate the psd of the white noise time-series data
3670 [a, S1, S2, navs] = prepare_analyze_noise_navs(noise_type, plist);
3671
3672 stest = true;
3673
3674 catch err
3675 disp(err.message)
3676 stest = false;
3677 end
3678 % </SyntaxCode>
3679
3680 % <AlgoDescription>
3681 %
3682 % 1) Check that calculated objects S1 and S2 are identical
3683 %
3684 % </AlgoDescription>
3685
3686 % <AlgoCode>
3687 atest = true;
3688
3689 if stest
3690 % Compare the output objects
3691 if ne(S1, S2, ple3)
3692 atest = false;
3693 end
3694 else
3695 atest = false;
3696 end
3697 % </AlgoCode>
3698
3699 % Return a result structure
3700 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3701 end % END UTP_52
3702
3703 %% UTP_53
3704
3705 % <TestDescription>
3706 %
3707 % Tests the possibility to set the number of averages rather than setting the Nfft:
3708 % 1) white noise produced from normal pdf, with:
3709 % a given mean value and sigma (distribution's 1st and 2nd order)
3710 % 2) psd of the noise, without detrending, random window, random navs
3711 % 3) get the number of averages
3712 % 4) get the nfft used
3713 % 5) run psd again, with the nfft used
3714 % 6) compare navs, nfft, psds
3715 %
3716
3717 % </TestDescription>
3718 function result = utp_53
3719
3720 % <SyntaxDescription>
3721 %
3722 % 1) white noise produced from normal pdf, with:
3723 % a given mean value and sigma (distribution's 1st and 2nd order)
3724 % 2) psd of the noise, without detrending, random window, random navs
3725 % 3) get the number of averages
3726 % 4) get the nfft used
3727 % 5) run psd again, with the nfft used
3728 % 6) run psd again, with conflicting parameters, and verify it uses
3729 % nfft rather than navs
3730 %
3731 % </SyntaxDescription>
3732
3733 % <SyntaxCode>
3734 try
3735 noise_type = 'Uniform';
3736
3737 % Evaluate the psd of the white noise time-series data
3738 [a, S1, S2, navs] = prepare_analyze_noise_navs(noise_type, plist);
3739
3740 npts_3 = fix(find(S1.hist.plistUsed, 'Nfft')/2);
3741
3742 % Calculates the psd asking for the number of points AND the window length
3743 pl_spec = S1.hist.plistUsed;
3744 pl_spec.pset('Nfft', npts_3, 'navs', navs);
3745 S3 = a.psd(pl_spec);
3746
3747 stest = true;
3748
3749 catch err
3750 disp(err.message)
3751 stest = false;
3752 end
3753 % </SyntaxCode>
3754
3755 % <AlgoDescription>
3756 %
3757 % 1) Check that calculated objects S1 and S2 are identical
3758 % 2) Check that S3 used different values
3759 %
3760 % </AlgoDescription>
3761
3762 % <AlgoCode>
3763 atest = true;
3764
3765 if stest
3766 % Compare the navs written in the output object with the requested one
3767 if ne(S1, S2, ple3) || ...
3768 ne(find(S3.hist.plistUsed, 'Nfft'), npts_3) || eq(S3.data.navs, navs)
3769 atest = false;
3770 end
3771 else
3772 atest = false;
3773 end
3774 % </AlgoCode>
3775
3776 % Return a result structure
3777 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3778 end % END UTP_53
3779
3780 %% UTP_54
3781
3782 % <TestDescription>
3783 %
3784 % Tests "conservation of energy":
3785 % 1) white noise produced from normal pdf, with a given mean value and
3786 % sigma (distribution's 1st and 2nd order)
3787 % 2) evaluate the sample mean value m and standard deviation s
3788 % 3) psd of the white noise
3789 % 4) compares the sqrt(sum(S * df)) with the standard deviation s
3790 %
3791
3792 % </TestDescription>
3793 function result = utp_54
3794
3795 % <SyntaxDescription>
3796 %
3797 % 1) Prepare the test tsdata:
3798 % white noise from normal distribution + offset
3799 % 2) Calculate the statistical parameters
3800 % 3) Estimate the psd
3801 %
3802 % </SyntaxDescription>
3803
3804 % <SyntaxCode>
3805 try
3806 % Array of parameters to pick from
3807 fs_list = [0.1;1;10];
3808 nsecs_list = [100:100:10000]';
3809 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
3810 mu_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
3811
3812 % Build time-series test data
3813
3814 % Picks the values at random from the list
3815 fs = utils.math.randelement(fs_list, 1);
3816 nsecs = utils.math.randelement(nsecs_list, 1);
3817 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
3818 mu_distr = utils.math.randelement(mu_distr_list, 1);
3819
3820 % White noise
3821 type = 'Normal';
3822 a_n = ao(plist('waveform', 'noise', ...
3823 'type', type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
3824 a_const = ao(mu_distr);
3825 a = a_n + a_const;
3826
3827 % Estimate the mean value and the sigma of the time-series data
3828 sigma_sample = std(a.y, 1);
3829 mu_sample = mean(a.y);
3830
3831 % Evaluate the psd of the time-series data
3832 % Parameters are crucial for the estimate to be correct!
3833 win = 'Rectangular';
3834 olap = 0;
3835 detrend = 0;
3836 n_pts = -1;
3837 scale = 'PSD';
3838
3839 S = a.psd(plist('Win', win, 'olap', olap, ...
3840 'Nfft', n_pts, 'order', detrend, 'scale', scale));
3841
3842 stest = true;
3843
3844 catch err
3845 disp(err.message)
3846 stest = false;
3847 end
3848 % </SyntaxCode>
3849
3850 % <AlgoDescription>
3851 %
3852 % 1) Check that calculated psd energy equals the rms content of the tsdata, estimated by the std of the sample
3853 %
3854 % </AlgoDescription>
3855
3856 % <AlgoCode>
3857 atest = true;
3858 TOL = 1e-12;
3859
3860 if stest
3861 % Integral of the spectrum
3862 sigma_psd = sqrt(sum(S.y * S.x(2)));
3863 if abs(sigma_psd - sigma_sample) / mean([sigma_psd sigma_sample]) >= TOL
3864 atest = false;
3865 end
3866 else
3867 atest = false;
3868 end
3869 % </AlgoCode>
3870
3871 % Return a result structure
3872 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3873 end % END UTP_54
3874
3875 %% UTP_60
3876
3877 % <TestDescription>
3878 %
3879 % Tests the 'basic' call:
3880 % 1) white noise produced from normal pdf, with:
3881 % a given mean value and sigma (distribution's 1st and 2nd order)
3882 % 2) psd of the noise, defualt detrending, default window, no averaging
3883 %
3884
3885 % </TestDescription>
3886 function result = utp_60
3887
3888 % <SyntaxDescription>
3889 %
3890 % 1) Prepare the test tsdata:
3891 % white noise from normal distribution + offset
3892 % 2) Calculate the statistical parameters
3893 % 3) Estimate the psd with default detrending, default window, no averaging
3894 %
3895 % </SyntaxDescription>
3896
3897 % <SyntaxCode>
3898 try
3899
3900 noise_type = 'Normal';
3901 win_type = char(prefs.getMiscPrefs.getDefaultWindow);
3902 scale = 'PSD';
3903
3904 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3905 plist('win_test', false, 'reduce_pts', false));
3906
3907 stest = true;
3908
3909 catch err
3910 disp(err.message)
3911 stest = false;
3912 end
3913 % </SyntaxCode>
3914
3915 % <AlgoDescription>
3916 %
3917 % 1) Nothing to check
3918 %
3919 % </AlgoDescription>
3920
3921 % <AlgoCode>
3922
3923 if stest
3924 atest = true;
3925 else
3926 atest = false;
3927 end
3928
3929 % </AlgoCode>
3930
3931 % Return a result structure
3932 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3933 end % END UTP_60
3934
3935 %% UTP_61
3936
3937 % <TestDescription>
3938 %
3939 % Tests the 'basic' call:
3940 % 1) white noise produced from normal pdf, with:
3941 % a given mean value and sigma (distribution's 1st and 2nd order)
3942 % 2) psd of the noise, defualt detrending, default window, no averaging
3943 %
3944
3945 % </TestDescription>
3946 function result = utp_61
3947
3948 % <SyntaxDescription>
3949 %
3950 % 1) Prepare the test tsdata:
3951 % white noise from normal distribution + offset
3952 % 2) Calculate the statistical parameters
3953 % 3) Estimate the psd with default detrending, default window, no averaging
3954 %
3955 % </SyntaxDescription>
3956
3957 % <SyntaxCode>
3958 try
3959
3960 noise_type = 'Normal';
3961 win_type = char(prefs.getMiscPrefs.getDefaultWindow);
3962 scale = 'ASD';
3963
3964 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
3965 plist('win_test', false, 'reduce_pts', false));
3966
3967 stest = true;
3968
3969 catch err
3970 disp(err.message)
3971 stest = false;
3972 end
3973 % </SyntaxCode>
3974
3975 % <AlgoDescription>
3976 %
3977 % 1) Nothing to check
3978 %
3979 % </AlgoDescription>
3980
3981 % <AlgoCode>
3982
3983 if stest
3984 atest = true;
3985 else
3986 atest = false;
3987 end
3988
3989 % </AlgoCode>
3990
3991 % Return a result structure
3992 result = utp_prepare_result(atest, stest, dbstack, mfilename);
3993 end % END UTP_61
3994
3995 %% UTP_62
3996
3997 % <TestDescription>
3998 %
3999 % Tests the 'basic' call:
4000 % 1) white noise produced from normal pdf, with:
4001 % a given mean value and sigma (distribution's 1st and 2nd order)
4002 % 2) psd of the noise, defualt detrending, default window, no averaging
4003 %
4004
4005 % </TestDescription>
4006 function result = utp_62
4007
4008 % <SyntaxDescription>
4009 %
4010 % 1) Prepare the test tsdata:
4011 % white noise from normal distribution + offset
4012 % 2) Calculate the statistical parameters
4013 % 3) Estimate the psd with default detrending, default window, no averaging
4014 %
4015 % </SyntaxDescription>
4016
4017 % <SyntaxCode>
4018 try
4019
4020 noise_type = 'Normal';
4021 win_type = char(prefs.getMiscPrefs.getDefaultWindow);
4022 scale = 'PS';
4023
4024 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
4025 plist('win_test', false, 'reduce_pts', false));
4026
4027 stest = true;
4028
4029 catch err
4030 disp(err.message)
4031 stest = false;
4032 end
4033 % </SyntaxCode>
4034
4035 % <AlgoDescription>
4036 %
4037 % 1) Nothing to check
4038 %
4039 % </AlgoDescription>
4040
4041 % <AlgoCode>
4042
4043 if stest
4044 atest = true;
4045 else
4046 atest = false;
4047 end
4048
4049 % </AlgoCode>
4050
4051 % Return a result structure
4052 result = utp_prepare_result(atest, stest, dbstack, mfilename);
4053 end % END UTP_62
4054
4055 %% UTP_63
4056
4057 % <TestDescription>
4058 %
4059 % Tests the 'basic' call:
4060 % 1) white noise produced from normal pdf, with:
4061 % a given mean value and sigma (distribution's 1st and 2nd order)
4062 % 2) psd of the noise, defualt detrending, default window, no averaging
4063 %
4064
4065 % </TestDescription>
4066 function result = utp_63
4067
4068 % <SyntaxDescription>
4069 %
4070 % 1) Prepare the test tsdata:
4071 % white noise from normal distribution + offset
4072 % 2) Calculate the statistical parameters
4073 % 3) Estimate the psd with default detrending, default window, no averaging
4074 %
4075 % </SyntaxDescription>
4076
4077 % <SyntaxCode>
4078 try
4079
4080 noise_type = 'Normal';
4081 win_type = char(prefs.getMiscPrefs.getDefaultWindow);
4082 scale = 'AS';
4083
4084 [a, S_name, S_obj, S_man] = prepare_analyze_noise_win(win_type, noise_type, scale, ...
4085 plist('win_test', false, 'reduce_pts', false));
4086
4087 stest = true;
4088
4089 catch err
4090 disp(err.message)
4091 stest = false;
4092 end
4093 % </SyntaxCode>
4094
4095 % <AlgoDescription>
4096 %
4097 % 1) Nothing to check
4098 %
4099 % </AlgoDescription>
4100
4101 % <AlgoCode>
4102
4103 if stest
4104 atest = true;
4105 else
4106 atest = false;
4107 end
4108
4109 % </AlgoCode>
4110
4111 % Return a result structure
4112 result = utp_prepare_result(atest, stest, dbstack, mfilename);
4113 end % END UTP_63
4114
4115 %% Helper function for window call construction
4116
4117 function [a, S_name, S_obj, S_w] = prepare_analyze_noise_win(win_type, noise_type, scale, pli)
4118 % Array of parameters to pick from
4119 fs_list = [0.1;1;2;5;10];
4120 nsecs_list = [20 100 1000:1000:10000]';
4121 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
4122 trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]';
4123
4124
4125 % Build time-series test data
4126
4127 % Picks the values at random from the list
4128 fs = utils.math.randelement(fs_list, 1);
4129 nsecs = utils.math.randelement(nsecs_list, 1);
4130 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
4131 trend_0 = utils.math.randelement(trend_0_list, 1);
4132
4133 % Pick units and prefix from those supported
4134 unit_list = unit.supportedUnits;
4135 % remove the first empty unit '' from the list, because then is it
4136 % possible that we add a prefix to an empty unit
4137 unit_list = unit_list(2:end);
4138 prefix_list = unit.supportedPrefixes;
4139
4140 % White noise
4141 a_n = ao(plist('waveform', 'noise', ...
4142 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
4143
4144 % Constant signal
4145 a_c = ao(trend_0);
4146
4147 % Total signal
4148 a = a_n + a_c;
4149
4150 % Set units
4151 a.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
4152
4153 % Evaluate the psd of the white noise time-series data
4154 n_pts_reduction_factor = find(pli, 'reduce_pts');
4155 if isempty(n_pts_reduction_factor) || ~n_pts_reduction_factor
4156 n_pts_reduction_factor = 1; % Taking single windows, no reduction in number of points
4157 end
4158 n_pts = a.nsecs * a.fs/ n_pts_reduction_factor;
4159
4160 olap = find(pli, 'olap');
4161 if isempty(olap)
4162 olap = 0; %% Should we take the rov instead?
4163 end
4164 if n_pts_reduction_factor == 1
4165 olap = 0; % This does not matter in the case of single window
4166 end
4167
4168 detrend_order = 0;
4169
4170 switch lower(win_type)
4171 case 'kaiser'
4172 psll = find(pli, 'psll');
4173 if isempty(psll)
4174 psll = find(ao.getInfo('psd').plists, 'psll');
4175 end
4176 win = specwin(win_type, fs*nsecs, psll);
4177 pl_spec = plist('Win', win_type, 'psll', psll, 'olap', olap, ...
4178 'Nfft', n_pts, 'order', detrend_order, 'scale', scale);
4179 otherwise
4180 win = specwin(win_type, fs*nsecs);
4181 pl_spec = plist('Win', win_type, 'olap', olap, ...
4182 'Nfft', n_pts, 'order', detrend_order, 'scale', scale);
4183 end
4184
4185 if find(pli, 'win_test')
4186 % Makes a copy of the data, and apply the window manually
4187 % Apply the detrending
4188 a_w = a.detrend(plist('N', detrend_order));
4189 % Apply the window, taking into account the correct normalization
4190 a_w.setY(a_w.y .* win.win'./sqrt(win.ws2 / win.len));
4191
4192 % Evaluate the psd of the windowed data
4193 S_w = a_w.psd(plist('Win', 'Rectangular', 'olap', olap, ...
4194 'Nfft', n_pts, 'order', -1, 'scale', scale));
4195
4196 % Calls psd applying the detrend and window internally
4197 % (passig window object)
4198 S_obj = a.psd(pl_spec);
4199
4200 else
4201 S_obj = ao;
4202 S_w = ao;
4203 end
4204
4205 % Calls psd applying the detrend and window internally
4206 % (passig window name)
4207 pl_spec.pset('Win', win_type);
4208 S_name = a.psd(pl_spec);
4209
4210 end
4211
4212
4213 %% Helper function for window call construction
4214
4215 function [sigma_sample, sigma_psd] = prepare_analyze_noise_energy(win_type, noise_type, scale, pli)
4216
4217 % Build time-series test data
4218 if find(pli, 'fs')
4219 fs = find(pli, 'fs');
4220 nsecs = find(pli, 'nsecs');
4221 sigma_distr = find(pli, 'sigma_distr');
4222 mu_distr = find(pli, 'mu_distr');
4223 else
4224 % Array of parameters to pick from
4225 fs_list = [0.1;1;10];
4226 nsecs_list = [100:100:10000];
4227 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
4228 mu_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
4229
4230 % Picks the values at random from the list
4231 fs = utils.math.randelement(fs_list, 1);
4232 nsecs = utils.math.randelement(nsecs_list, 1);
4233 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
4234 mu_distr = utils.math.randelement(mu_distr_list, 1);
4235 end
4236
4237 % Pick units and prefix from those supported
4238 unit_list = unit.supportedUnits;
4239 % remove the first empty unit '' from the list, because then is it
4240 % possible that we add a prefix to an empty unit
4241 unit_list = unit_list(2:end);
4242 prefix_list = unit.supportedPrefixes;
4243
4244 % White noise
4245 a_n = ao(plist('waveform', 'noise', ...
4246 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
4247
4248 % Constant signal
4249 a_c = ao(mu_distr);
4250
4251 % Total signal
4252 a = a_n + a_c;
4253
4254 % Set units
4255 a.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
4256
4257 % Evaluate the psd of the white noise time-series data
4258 olap = find(pli, 'olap');
4259 if isempty(olap)
4260 olap = 0; %% Should we take the rov instead?
4261 end
4262
4263 detrend_order = find(pli, 'detrend_order');
4264 if isempty(detrend_order)
4265 detrend_order = 0;
4266 end
4267
4268 n_pts = -1;
4269
4270 switch lower(win_type)
4271 case 'kaiser'
4272 psll = find(pli, 'psll');
4273 if isempty(psll)
4274 psll = find(ao.getInfo('psd').plists, 'psll');
4275 end
4276 win = specwin(win_type, fs*nsecs, psll);
4277 pl_spec = plist('Win', win_type, 'psll', psll, 'olap', olap, ...
4278 'Nfft', n_pts, 'order', detrend_order, 'scale', scale);
4279 otherwise
4280 win = specwin(win_type, fs*nsecs);
4281 pl_spec = plist('Win', win_type, 'olap', olap, ...
4282 'Nfft', n_pts, 'order', detrend_order, 'scale', scale);
4283 end
4284
4285 % Calls psd applying the detrend and window internally
4286 % (passig window name)
4287 pl_spec.pset('Win', win_type);
4288 S = a.psd(pl_spec);
4289
4290 % Evaluate the statistics of the time-series and spectra
4291
4292 % Estimate the mean value and the sigma of the time-series data
4293 sigma_sample = std(a.y, 1);
4294 mu_sample = mean(a.y);
4295
4296 % Integral of the spectrum
4297 sigma_psd = sqrt(sum(S.y * S.x(2)));
4298
4299 end
4300
4301 %% Helper function for window call construction, navs option
4302
4303 function [a, S, S1, navs] = prepare_analyze_noise_navs(noise_type, pli)
4304 % Array of parameters to pick from
4305 fs_list = [0.1;1;2;5;10];
4306 nsecs_list = [2000:1000:10000]';
4307 sigma_distr_list = [1e-6 2e-3 0.25 1:0.1:10]';
4308 trend_0_list = [1e-6 2e-3 0.25 1:0.1:10]';
4309
4310
4311 % Build time-series test data
4312
4313 % Picks the values at random from the list
4314 fs = utils.math.randelement(fs_list, 1);
4315 nsecs = utils.math.randelement(nsecs_list, 1);
4316 sigma_distr = utils.math.randelement(sigma_distr_list, 1);
4317 trend_0 = utils.math.randelement(trend_0_list, 1);
4318
4319 % Pick units and prefix from those supported
4320 unit_list = unit.supportedUnits;
4321 % remove the first empty unit '' from the list, because then is it
4322 % possible that we add a prefix to an empty unit
4323 unit_list = unit_list(2:end);
4324 prefix_list = unit.supportedPrefixes;
4325
4326 % White noise
4327 a_n = ao(plist('waveform', 'noise', ...
4328 'type', noise_type, 'fs', fs, 'nsecs', nsecs, 'sigma', sigma_distr));
4329
4330 % Constant signal
4331 a_c = ao(trend_0);
4332
4333 % Total signals
4334 a = a_n + a_c;
4335
4336 % Set units
4337 a.setYunits(unit([cell2mat(utils.math.randelement(prefix_list,1)) cell2mat(utils.math.randelement(unit_list,1))]));
4338
4339 % Evaluate the psd of the white noise time-series data
4340 olap = 0;
4341 detrend_order = 0;
4342 n_pts = -1;
4343
4344 navs = fix(utils.math.randelement(logspace(0,log10(max(0,a.len/10)),50),1));
4345
4346 % Evaluate the psd of the white noise time-series data
4347 % Window
4348 win_list = specwin.getTypes;
4349 win_type = utils.math.randelement(win_list,1);
4350 win_type = win_type{1};
4351
4352 switch lower(win_type)
4353 case 'kaiser'
4354 psll = utils.math.randelement([0:10:200],1);
4355 if psll == 0
4356 psll = find(ao.getInfo('psd').plists, 'psll');
4357 end
4358 pl_spec = plist('Win', win_type, 'psll', psll, 'olap', olap, 'order', detrend_order);
4359 case 'levelledhanning'
4360 pl_spec = plist('Win', 'BH92', 'olap', olap, 'order', detrend_order);
4361 otherwise
4362 pl_spec = plist('Win', win_type, 'olap', olap, 'order', detrend_order);
4363 end
4364
4365 % Calls psd asking for the number of averages
4366 pl_spec.pset('Nfft', n_pts, 'navs', navs);
4367 S = a.psd(pl_spec);
4368
4369 % Calls psd asking for the number of points just evaluated
4370 pl_spec.pset('Nfft', find(S.hist.plistUsed, 'Nfft'));
4371 pl_spec.remove('navs');
4372 S1 = a.psd(pl_spec);
4373
4374 end
4375
4376 %% Helper function for window call construction, navs option
4377 function [S_ap, S_sigma_total, S_expected_total] = test_psd_energy_prepare_data(win_type, detrend_order, pli)
4378
4379 % Calculate the expected level of noise from sample statistics
4380 fs = a_n_long.fs;
4381 S_sigma_total = std(a_n_long.y, 1);
4382 S_expected_total = S_sigma_total^2 / fs * 2;
4383
4384 olap = 0;
4385 n_pts = -1;
4386 scale = 'PSD';
4387
4388 switch lower(win_type)
4389 case 'kaiser'
4390 psll = find(pli, 'psll');
4391 if isempty(psll)
4392 psll = find(ao.getInfo('psd').plists, 'psll');
4393 end
4394 pl_spec = plist('Win', win_type, 'psll', psll, 'olap', olap, ...
4395 'Nfft', n_pts, 'order', detrend_order, 'scale', scale);
4396 case 'levelledhanning'
4397 pl_spec = plist('Win', 'BH92', 'olap', olap, ...
4398 'Nfft', n_pts, 'order', detrend_order, 'scale', scale);
4399 otherwise
4400 pl_spec = plist('Win', win_type, 'olap', olap, ...
4401 'Nfft', n_pts, 'order', detrend_order, 'scale', scale);
4402 end
4403
4404 % Calls psd applying the detrend and window internally
4405 % (passig window name)
4406 S_ap = ap.psd(pl_spec);
4407
4408 end
4409
4410 %% Algorithm test for window comparison
4411
4412 function varargout = algo_test_psd_win(S_name, S_obj, S_man, pli)
4413 atest = true;
4414 msg = '';
4415 pl_def = plist('TOL', 1e-12);
4416 TOL = find(parse(pli, pl_def), 'TOL');
4417
4418 % Compare the psd evaluated with the window name / object
4419 if ~eq(S_name, S_obj, ple1)
4420 atest = false;
4421 msg = ['Two PSD objects are not equal: ' lastwarn];
4422 end
4423 % Compare the psd evaluated with the two methods
4424 res = abs(S_obj - S_man) ./ S_obj;
4425 if any(res.y >= TOL)
4426 atest = false;
4427 msg = 'Two PSDs are not equal';
4428 end
4429
4430 if nargout == 1
4431 varargout{1} = atest;
4432 elseif nargout == 2
4433 varargout{1} = atest;
4434 varargout{2} = msg;
4435 else
4436 error('Incorrect outputs');
4437 end
4438
4439 end
4440
4441
4442 %% Algorithm test for energy conservation tests
4443
4444 function atest = algo_test_psd_energy(ap, S_ap, S_expected_total, pli)
4445
4446 % Checks on individual PSDs: mean and standard deviation of the
4447 % PSD points
4448 fs = ap.fs;
4449 Npoints = length(S_ap(1).y);
4450 S_mean = mean(S_ap.y);
4451 S_err = std(S_ap.y, 1) / sqrt(Npoints);
4452
4453 % Expected value, estimated from the std of the individual segments
4454 S_sigma = std(ap.y, 1);
4455 S_expected = S_sigma.^2 / fs * 2;
4456
4457 % Comparison with the mean performed on the individual segments
4458 test_level = abs(S_mean - S_expected) < S_err;
4459 sum(test_level) / Nsegments; % TODO: implement hypothesis test here
4460
4461 % Checks on averaged PSDs: mean and standard deviation of all the
4462 % PSD points
4463 S_mean_total = mean(S_mean);
4464 S_err_total = std(S_mean) / sqrt(Nsegments);
4465
4466 % Compares the grand-mean with the estimated white noise level
4467 if sum(abs(S_mean_total - S_expected_total) < S_err_total)
4468 atest = true;
4469 else
4470 atest = false;
4471 end
4472
4473 end
4474
4475
4476 end