Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/psdconf.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % PSDCONF Calculates confidence levels and variance for psd
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: PSDCONF Input a spectrum estimated with psd or lpsd (Welch's
5 % Overlapped Segment Averaging Method) and calculates confidence levels and
6 % variance for the spectrum.
7 %
8 % CALL: [lcl,ucl] = psdconf(a,pl)
9 % [lcl,ucl,var] = psdconf(a,pl)
10 %
11 % INPUTS:
12 % a - input analysis objects containing power spectral
13 % densities calculated with psd or lpsd.
14 % pl - input parameter list
15 %
16 % OUTPUTS:
17 % lcl - lower confidence level
18 % ucl - upper confidence level
19 % var - expected spectrum variance
20 %
21 %
22 % If the last input argument is a parameter list (plist).
23 % The following parameters are recognised.
24 %
25 %
26 % NOTE1: PSDCONF checks the navs field to distinguish between psd and lpsd
27 % power spectra. If a.data.navs is NaN then it assumes to dealing with lpsd
28 % power spectrum at input.
29 %
30 % NOTE2: Copied directly from MATLAB chi2conf function (Signal Processing
31 % Toolbox) and extended to do degrees of fredom calculation, variance
32 % calculation, to input AOs and to input plist.
33 % Copyright 2007 The MathWorks, Inc.
34 % Revision: 1.6.4.4 Date: 2008/05/31 23:27:28
35 %
36 %
37 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'psdconf')">Parameters Description</a>
38 %
39 % VERSION: $Id: psdconf.m,v 1.14 2011/05/18 07:53:57 mauro Exp $
40 %
41 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42
43 function varargout = psdconf(varargin)
44
45 %%% Check if this is a call for parameters
46 if utils.helper.isinfocall(varargin{:})
47 varargout{1} = getInfo(varargin{3});
48 return
49 end
50
51 import utils.const.*
52 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
53
54 %%% Collect input variable names
55 in_names = cell(size(varargin));
56 for ii = 1:nargin,in_names{ii} = inputname(ii);end
57
58 %%% Collect all AOs
59 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
60 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
61
62 %%% avoid multiple AO at input
63 if numel(as)>1
64 error('!!! Too many input AOs, PSDCONF can process only one AO per time !!!')
65 end
66
67 %%% check that fsdata is input
68 if ~isa(as.data, 'fsdata')
69 error('!!! Non-fsdata input, PSDCONF can process only fsdata !!!')
70 end
71
72 %%% avoid input modification
73 if nargout == 0
74 error('!!! PSDCONF cannot be used as a modifier. Please give an output variable !!!');
75 end
76
77 %%% Parse plists
78 pl = parse(pl, getDefaultPlist());
79
80 %%% Find parameters
81 conf = find(pl, 'conf');
82 conf = conf/100; % go from percentage to fractional
83
84 %%% getting metho
85 if isnan(as.data.navs)
86 mtd = 'lpsd';
87 else
88 mtd = 'psd';
89 end
90
91 %%% switching over methods
92 switch mtd
93 case 'psd'
94 %%% confidence levels for spectra calculated with psd
95 % get number of averages
96 navs = as.data.navs;
97 % get window object
98 w = as.hist.plistUsed.find('WIN');
99 % percentage of overlap
100 olap = as.hist.plistUsed.find('OLAP')./100;
101 % number of bins in each fft
102 nfft = as.hist.plistUsed.find('NFFT');
103
104 % Normalize window data in order to be square integrable to 1
105 win = w.win ./ sqrt(w.ws2);
106
107 % Calculates total number of data in the original time-series
108 Ntot = ceil(navs*(nfft-olap*nfft)+olap*nfft);
109
110 % defining the shift factor
111 n = floor((Ntot-nfft)./(navs-1));
112
113 % Calculating correction factor
114 sm = 0;
115 % Padding to zero the window
116 % willing to work with columns
117 [ii,kk] = size(win);
118 if ii<kk
119 win = win.';
120 end
121 win = [win; zeros(navs*n,1)];
122 for mm = 1:navs-1
123 pc = 0;
124 for tt = 1:nfft
125 pc = pc + win(tt)*win(tt + mm*n);
126 end
127 sm = sm + (1 - mm/navs)*(abs(pc)^2);
128 end
129 % Calculating degrees of freedom
130 dof = (2*navs)/(1+2*sm);
131 dof = round(dof);
132
133 % Calculating Confidence Levels factors
134 alfa = 1 - conf;
135 c = chi2inv([1-alfa/2 alfa/2],dof);
136 c = dof./c;
137
138 % calculating variance
139 expvar = ((as.data.y).^2).*2./dof;
140
141 % calculating confidence levels
142 lwb = as.data.y.*c(1);
143 upb = as.data.y.*c(2);
144
145 case 'lpsd'
146 %%% confidence levels for spectra calculated with lpsd
147 % get window used
148 uwin = as.hist.plistUsed.find('WIN');
149
150 % extract number of frequencies bins
151 nf = length(as.x);
152
153 % dft length for each bin
154 if ~isempty(as.procinfo)
155 L = as.procinfo.find('L');
156 else
157 error('### The AO doesn''t have any procinfo with the key ''L''');
158 end
159
160 % set original data length as the length of the first window
161 nx = L(1);
162
163 % windows overlap
164 olap = as.hist.plistUsed.find('OLAP');
165
166 dofs = ones(nf,1);
167 cl = ones(nf,2);
168 for jj=1:nf
169 l = L(jj);
170 % compute window
171 switch uwin.type
172 case 'Kaiser'
173 w = specwin(uwin.type, l, uwin.psll);
174 otherwise
175 w = specwin(uwin.type, l);
176 end
177 % Normalize window data in order to be square integrable to 1
178 owin = w.win ./ sqrt(w.ws2);
179
180 % Compute the number of averages we want here
181 segLen = l;
182 nData = nx;
183 ovfact = 1 / (1 - olap);
184 davg = (((nData - segLen)) * ovfact) / segLen + 1;
185 navg = round(davg);
186
187 % Compute steps between segments
188 if navg == 1
189 shift = 1;
190 else
191 shift = (nData - segLen) / (navg - 1);
192 end
193 if shift < 1 || isnan(shift)
194 shift = 1;
195 end
196 n = floor(shift);
197
198 % Calculating correction factor
199 sm = 0;
200 % Padding to zero the window
201 % willing to work with columns
202 [ii,kk] = size(owin);
203 if ii<kk
204 owin = owin.';
205 end
206 win = [owin; zeros(navg*n,1)];
207 for mm = 1:navg-1
208 pc = 0;
209 for tt = 1:segLen
210 pc = pc + win(tt)*win(tt + mm*n);
211 end
212 sm = sm + (1 - mm/navg)*(abs(pc)^2);
213 end
214 % Calculating degrees of freedom
215 dof = (2*navg)/(1+2*sm);
216 dof = round(dof);
217
218 % Calculating Confidence Levels factors
219 alfa = 1 - conf;
220 c = chi2inv([1-alfa/2 alfa/2],dof);
221 c = dof./c;
222
223 % storing c and dof
224 dofs(jj) = dof;
225 cl(jj,1) = c(1);
226 cl(jj,2) = c(2);
227 end % for jj=1:nf
228
229 % willing to work with columns
230 dy = as.data.y;
231 [ii,kk] = size(dy);
232 if ii<kk
233 dy = dy.';
234 rsp = true;
235 else
236 rsp = false;
237 end
238 % calculating variance
239 expvar = ((dy).^2).*2./dofs;
240
241 % calculating confidence levels
242 lwb = dy.*cl(:,1);
243 upb = dy.*cl(:,2);
244
245 % reshaping if necessary
246 if rsp
247 expvar = expvar.';
248 lwb = lwb.';
249 upb = upb.';
250 end
251
252 end %switch mtd
253
254 % Output data
255
256 % defining units
257 inputunit = get(as.data,'yunits');
258 varunit = unit(inputunit.^2);
259 varunit.simplify;
260 levunit = inputunit;
261 levunit.simplify;
262
263
264 % variance
265 dvar = fsdata();
266 dvar.setFs(as.data.fs);
267 dvar.setT0(as.data.t0);
268 dvar.setEnbw(as.data.enbw);
269 dvar.setNavs(as.data.navs);
270 dvar.setXunits(copy(as.data.xunits,1));
271 dvar.setYunits(varunit);
272 dvar.setX(as.data.x);
273 dvar.setY(expvar);
274 ovar = ao(dvar);
275 % Set output AO name
276 ovar.name = sprintf('var(%s)', ao_invars{:});
277 % Add history
278 ovar.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
279
280 % lower confidence level
281 dlwb = fsdata();
282 dlwb.setFs(as.data.fs);
283 dlwb.setT0(as.data.t0);
284 dlwb.setEnbw(as.data.enbw);
285 dlwb.setNavs(as.data.navs);
286 dlwb.setXunits(copy(as.data.xunits,1));
287 dlwb.setYunits(levunit);
288 dlwb.setX(as.data.x);
289 dlwb.setY(lwb);
290 olwb = ao(dlwb);
291 % Set output AO name
292 clev = [num2str(conf*100) '%'];
293 olwb.name = sprintf('%s_low_conf_level(%s)', clev, ao_invars{:});
294 % Add history
295 olwb.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
296
297 % upper confidence level
298 dupb = fsdata();
299 dupb.setFs(as.data.fs);
300 dupb.setT0(as.data.t0);
301 dupb.setEnbw(as.data.enbw);
302 dupb.setNavs(as.data.navs);
303 dupb.setXunits(copy(as.data.xunits,1));
304 dupb.setYunits(levunit);
305 dupb.setX(as.data.x);
306 dupb.setY(upb);
307 oupb = ao(dupb);
308 % Set output AO name
309 oupb.name = sprintf('%s_up_conf_level(%s)', clev, ao_invars{:});
310 % Add history
311 oupb.addHistory(getInfo('None'), pl, [ao_invars(:)], [as.hist]);
312
313 % output
314 if nargout == 2
315 varargout{1} = olwb; % lower conf level
316 varargout{2} = oupb; % upper conf level
317 elseif nargout == 3
318 varargout{1} = olwb; % lower conf level
319 varargout{2} = oupb; % upper conf level
320 varargout{3} = ovar; % expected variance
321 end
322
323 end
324
325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
326 % Local Functions %
327 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
328
329 %----- Matlab Functions ---------------------------------------------------
330
331 %--------------------------------------------------------------------------
332 function x = chi2inv(p,v)
333 %CHI2INV Inverse of the chi-square cumulative distribution function (cdf).
334 % X = CHI2INV(P,V) returns the inverse of the chi-square cdf with V
335 % degrees of freedom at the values in P. The chi-square cdf with V
336 % degrees of freedom, is the gamma cdf with parameters V/2 and 2.
337 %
338 % The size of X is the common size of P and V. A scalar input
339 % functions as a constant matrix of the same size as the other input.
340
341 % References:
342 % [1] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
343 % Functions", Government Printing Office, 1964, 26.4.
344 % [2] E. Kreyszig, "Introductory Mathematical Statistics",
345 % John Wiley, 1970, section 10.2 (page 144)
346
347
348 if nargin < 2,
349 error(generatemsgid('Nargchk'),'Requires two input arguments.');
350 end
351
352 [errorcode p v] = distchck(2,p,v);
353
354 if errorcode > 0
355 error(generatemsgid('InvalidDimensions'),'Requires non-scalar arguments to match in size.');
356 end
357
358 % Call the gamma inverse function.
359 x = gaminv(p,v/2,2);
360
361 % Return NaN if the degrees of freedom is not a positive integer.
362 k = find(v < 0 | round(v) ~= v);
363 if any(k)
364 tmp = NaN;
365 x(k) = tmp(ones(size(k)));
366 end
367
368 end
369
370 %--------------------------------------------------------------------------
371
372 function x = gaminv(p,a,b)
373 %GAMINV Inverse of the gamma cumulative distribution function (cdf).
374 % X = GAMINV(P,A,B) returns the inverse of the gamma cdf with
375 % parameters A and B, at the probabilities in P.
376 %
377 % The size of X is the common size of the input arguments. A scalar input
378 % functions as a constant matrix of the same size as the other inputs.
379 %
380 % GAMINV uses Newton's method to converge to the solution.
381
382 % References:
383 % [1] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
384 % Functions", Government Printing Office, 1964, 6.5.
385
386 % B.A. Jones 1-12-93
387 % Was: Revision: 1.2, Date: 1996/07/25 16:23:36
388
389 if nargin<3,
390 b=1;
391 end
392
393 [errorcode p a b] = distchck(3,p,a,b);
394
395 if errorcode > 0
396 error(generatemsgid('InvalidDimensions'),'The arguments must be the same size or be scalars.');
397 end
398
399 % Initialize X to zero.
400 x = zeros(size(p));
401
402 k = find(p<0 | p>1 | a <= 0 | b <= 0);
403 if any(k),
404 tmp = NaN;
405 x(k) = tmp(ones(size(k)));
406 end
407
408 % The inverse cdf of 0 is 0, and the inverse cdf of 1 is 1.
409 k0 = find(p == 0 & a > 0 & b > 0);
410 if any(k0),
411 x(k0) = zeros(size(k0));
412 end
413
414 k1 = find(p == 1 & a > 0 & b > 0);
415 if any(k1),
416 tmp = Inf;
417 x(k1) = tmp(ones(size(k1)));
418 end
419
420 % Newton's Method
421 % Permit no more than count_limit iterations.
422 count_limit = 100;
423 count = 0;
424
425 k = find(p > 0 & p < 1 & a > 0 & b > 0);
426 pk = p(k);
427
428 % Supply a starting guess for the iteration.
429 % Use a method of moments fit to the lognormal distribution.
430 mn = a(k) .* b(k);
431 v = mn .* b(k);
432 temp = log(v + mn .^ 2);
433 mu = 2 * log(mn) - 0.5 * temp;
434 sigma = -2 * log(mn) + temp;
435 xk = exp(norminv(pk,mu,sigma));
436
437 h = ones(size(pk));
438
439 % Break out of the iteration loop for three reasons:
440 % 1) the last update is very small (compared to x)
441 % 2) the last update is very small (compared to sqrt(eps))
442 % 3) There are more than 100 iterations. This should NEVER happen.
443
444 while(any(abs(h) > sqrt(eps)*abs(xk)) & max(abs(h)) > sqrt(eps) ...
445 & count < count_limit),
446
447 count = count + 1;
448 h = (gamcdf(xk,a(k),b(k)) - pk) ./ gampdf(xk,a(k),b(k));
449 xnew = xk - h;
450 % Make sure that the current guess stays greater than zero.
451 % When Newton's Method suggests steps that lead to negative guesses
452 % take a step 9/10ths of the way to zero:
453 ksmall = find(xnew < 0);
454 if any(ksmall),
455 xnew(ksmall) = xk(ksmall) / 10;
456 h = xk-xnew;
457 end
458 xk = xnew;
459 end
460
461
462 % Store the converged value in the correct place
463 x(k) = xk;
464
465 if count == count_limit,
466 fprintf('\nWarning: GAMINV did not converge.\n');
467 str = 'The last step was: ';
468 outstr = sprintf([str,'%13.8f'],h);
469 fprintf(outstr);
470 end
471 end
472
473 %--------------------------------------------------------------------------
474 function x = norminv(p,mu,sigma)
475 %NORMINV Inverse of the normal cumulative distribution function (cdf).
476 % X = NORMINV(P,MU,SIGMA) finds the inverse of the normal cdf with
477 % mean, MU, and standard deviation, SIGMA.
478 %
479 % The size of X is the common size of the input arguments. A scalar input
480 % functions as a constant matrix of the same size as the other inputs.
481 %
482 % Default values for MU and SIGMA are 0 and 1 respectively.
483
484 % References:
485 % [1] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
486 % Functions", Government Printing Office, 1964, 7.1.1 and 26.2.2
487
488 % Was: Revision: 1.2, Date: 1996/07/25 16:23:36
489
490 if nargin < 3,
491 sigma = 1;
492 end
493
494 if nargin < 2;
495 mu = 0;
496 end
497
498 [errorcode p mu sigma] = distchck(3,p,mu,sigma);
499
500 if errorcode > 0
501 error(generatemsgid('InvalidDimensions'),'Requires non-scalar arguments to match in size.');
502 end
503
504 % Allocate space for x.
505 x = zeros(size(p));
506
507 % Return NaN if the arguments are outside their respective limits.
508 k = find(sigma <= 0 | p < 0 | p > 1);
509 if any(k)
510 tmp = NaN;
511 x(k) = tmp(ones(size(k)));
512 end
513
514 % Put in the correct values when P is either 0 or 1.
515 k = find(p == 0);
516 if any(k)
517 tmp = Inf;
518 x(k) = -tmp(ones(size(k)));
519 end
520
521 k = find(p == 1);
522 if any(k)
523 tmp = Inf;
524 x(k) = tmp(ones(size(k)));
525 end
526
527 % Compute the inverse function for the intermediate values.
528 k = find(p > 0 & p < 1 & sigma > 0);
529 if any(k),
530 x(k) = sqrt(2) * sigma(k) .* erfinv(2 * p(k) - 1) + mu(k);
531 end
532 end
533
534 %--------------------------------------------------------------------------
535 function p = gamcdf(x,a,b)
536 %GAMCDF Gamma cumulative distribution function.
537 % P = GAMCDF(X,A,B) returns the gamma cumulative distribution
538 % function with parameters A and B at the values in X.
539 %
540 % The size of P is the common size of the input arguments. A scalar input
541 % functions as a constant matrix of the same size as the other inputs.
542 %
543 % Some references refer to the gamma distribution with a single
544 % parameter. This corresponds to the default of B = 1.
545 %
546 % GAMMAINC does computational work.
547
548 % References:
549 % [1] L. Devroye, "Non-Uniform Random Variate Generation",
550 % Springer-Verlag, 1986. p. 401.
551 % [2] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
552 % Functions", Government Printing Office, 1964, 26.1.32.
553
554 % Was: Revision: 1.2, Date: 1996/07/25 16:23:36
555
556 if nargin < 3,
557 b = 1;
558 end
559
560 if nargin < 2,
561 error(generatemsgid('Nargchk'),'Requires at least two input arguments.');
562 end
563
564 [errorcode x a b] = distchck(3,x,a,b);
565
566 if errorcode > 0
567 error(generatemsgid('InvalidDimensions'),'Requires non-scalar arguments to match in size.');
568 end
569
570 % Return NaN if the arguments are outside their respective limits.
571 k1 = find(a <= 0 | b <= 0);
572 if any(k1)
573 tmp = NaN;
574 p(k1) = tmp(ones(size(k1)));
575 end
576
577 % Initialize P to zero.
578 p = zeros(size(x));
579
580 k = find(x > 0 & ~(a <= 0 | b <= 0));
581 if any(k),
582 p(k) = gammainc(x(k) ./ b(k),a(k));
583 end
584
585 % Make sure that round-off errors never make P greater than 1.
586 k = find(p > 1);
587 if any(k)
588 p(k) = ones(size(k));
589 end
590 end
591
592 %--------------------------------------------------------------------------
593 function y = gampdf(x,a,b)
594 %GAMPDF Gamma probability density function.
595 % Y = GAMPDF(X,A,B) returns the gamma probability density function
596 % with parameters A and B, at the values in X.
597 %
598 % The size of Y is the common size of the input arguments. A scalar input
599 % functions as a constant matrix of the same size as the other inputs.
600 %
601 % Some references refer to the gamma distribution with a single
602 % parameter. This corresponds to the default of B = 1.
603
604 % References:
605 % [1] L. Devroye, "Non-Uniform Random Variate Generation",
606 % Springer-Verlag, 1986, pages 401-402.
607
608 % Was: Revision: 1.2, Date: 1996/07/25 16:23:36
609
610 if nargin < 3,
611 b = 1;
612 end
613
614 if nargin < 2,
615 error(generatemsgid('Nargchk'),'Requires at least two input arguments');
616 end
617
618 [errorcode x a b] = distchck(3,x,a,b);
619
620 if errorcode > 0
621 error(generatemsgid('InvalidDimensions'),'Requires non-scalar arguments to match in size.');
622 end
623
624 % Initialize Y to zero.
625 y = zeros(size(x));
626
627 % Return NaN if the arguments are outside their respective limits.
628 k1 = find(a <= 0 | b <= 0);
629 if any(k1)
630 tmp = NaN;
631 y(k1) = tmp(ones(size(k1)));
632 end
633
634 k=find(x > 0 & ~(a <= 0 | b <= 0));
635 if any(k)
636 y(k) = (a(k) - 1) .* log(x(k)) - (x(k) ./ b(k)) - gammaln(a(k)) - a(k) .* log(b(k));
637 y(k) = exp(y(k));
638 end
639 k1 = find(x == 0 & a < 1);
640 if any(k1)
641 tmp = Inf;
642 y(k1) = tmp(ones(size(k1)));
643 end
644 k2 = find(x == 0 & a == 1);
645 if any(k2)
646 y(k2) = (1./b(k2));
647 end
648 end
649
650 function [errorcode,out1,out2,out3,out4] = distchck(nparms,arg1,arg2,arg3,arg4)
651 %DISTCHCK Checks the argument list for the probability functions.
652
653 % B.A. Jones 1-22-93
654 % Was: Revision: 1.2, Date: 1996/07/25 16:23:36
655
656 errorcode = 0;
657
658 if nparms == 1
659 out1 = arg1;
660 return;
661 end
662
663 if nparms == 2
664 [r1 c1] = size(arg1);
665 [r2 c2] = size(arg2);
666 scalararg1 = (prod(size(arg1)) == 1);
667 scalararg2 = (prod(size(arg2)) == 1);
668 if ~scalararg1 & ~scalararg2
669 if r1 ~= r2 | c1 ~= c2
670 errorcode = 1;
671 return;
672 end
673 end
674 if scalararg1
675 out1 = arg1(ones(r2,1),ones(c2,1));
676 else
677 out1 = arg1;
678 end
679 if scalararg2
680 out2 = arg2(ones(r1,1),ones(c1,1));
681 else
682 out2 = arg2;
683 end
684 end
685
686 if nparms == 3
687 [r1 c1] = size(arg1);
688 [r2 c2] = size(arg2);
689 [r3 c3] = size(arg3);
690 scalararg1 = (prod(size(arg1)) == 1);
691 scalararg2 = (prod(size(arg2)) == 1);
692 scalararg3 = (prod(size(arg3)) == 1);
693
694 if ~scalararg1 & ~scalararg2
695 if r1 ~= r2 | c1 ~= c2
696 errorcode = 1;
697 return;
698 end
699 end
700
701 if ~scalararg1 & ~scalararg3
702 if r1 ~= r3 | c1 ~= c3
703 errorcode = 1;
704 return;
705 end
706 end
707
708 if ~scalararg3 & ~scalararg2
709 if r3 ~= r2 | c3 ~= c2
710 errorcode = 1;
711 return;
712 end
713 end
714
715 if ~scalararg1
716 out1 = arg1;
717 end
718 if ~scalararg2
719 out2 = arg2;
720 end
721 if ~scalararg3
722 out3 = arg3;
723 end
724 rows = max([r1 r2 r3]);
725 columns = max([c1 c2 c3]);
726
727 if scalararg1
728 out1 = arg1(ones(rows,1),ones(columns,1));
729 end
730 if scalararg2
731 out2 = arg2(ones(rows,1),ones(columns,1));
732 end
733 if scalararg3
734 out3 = arg3(ones(rows,1),ones(columns,1));
735 end
736 out4 =[];
737
738 end
739
740 if nparms == 4
741 [r1 c1] = size(arg1);
742 [r2 c2] = size(arg2);
743 [r3 c3] = size(arg3);
744 [r4 c4] = size(arg4);
745 scalararg1 = (prod(size(arg1)) == 1);
746 scalararg2 = (prod(size(arg2)) == 1);
747 scalararg3 = (prod(size(arg3)) == 1);
748 scalararg4 = (prod(size(arg4)) == 1);
749
750 if ~scalararg1 & ~scalararg2
751 if r1 ~= r2 | c1 ~= c2
752 errorcode = 1;
753 return;
754 end
755 end
756
757 if ~scalararg1 & ~scalararg3
758 if r1 ~= r3 | c1 ~= c3
759 errorcode = 1;
760 return;
761 end
762 end
763
764 if ~scalararg1 & ~scalararg4
765 if r1 ~= r4 | c1 ~= c4
766 errorcode = 1;
767 return;
768 end
769 end
770
771 if ~scalararg3 & ~scalararg2
772 if r3 ~= r2 | c3 ~= c2
773 errorcode = 1;
774 return;
775 end
776 end
777
778 if ~scalararg4 & ~scalararg2
779 if r4 ~= r2 | c4 ~= c2
780 errorcode = 1;
781 return;
782 end
783 end
784
785 if ~scalararg3 & ~scalararg4
786 if r3 ~= r4 | c3 ~= c4
787 errorcode = 1;
788 return;
789 end
790 end
791
792
793 if ~scalararg1
794 out1 = arg1;
795 end
796 if ~scalararg2
797 out2 = arg2;
798 end
799 if ~scalararg3
800 out3 = arg3;
801 end
802 if ~scalararg4
803 out4 = arg4;
804 end
805
806 rows = max([r1 r2 r3 r4]);
807 columns = max([c1 c2 c3 c4]);
808 if scalararg1
809 out1 = arg1(ones(rows,1),ones(columns,1));
810 end
811 if scalararg2
812 out2 = arg2(ones(rows,1),ones(columns,1));
813 end
814 if scalararg3
815 out3 = arg3(ones(rows,1),ones(columns,1));
816 end
817 if scalararg4
818 out4 = arg4(ones(rows,1),ones(columns,1));
819 end
820 end
821 end
822
823 %----- LTPDA FUNCTIONS ----------------------------------------------------
824 %--------------------------------------------------------------------------
825 % Get Info Object
826 %--------------------------------------------------------------------------
827 function ii = getInfo(varargin)
828 if nargin == 1 && strcmpi(varargin{1}, 'None')
829 sets = {};
830 pl = [];
831 else
832 sets = {'Default'};
833 pl = getDefaultPlist;
834 end
835 % Build info object
836 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: psdconf.m,v 1.14 2011/05/18 07:53:57 mauro Exp $', sets, pl);
837 ii.setModifier(false);
838 ii.setOutmin(2);
839 ii.setOutmax(3);
840 end
841
842 %--------------------------------------------------------------------------
843 % Get Default Plist
844 %--------------------------------------------------------------------------
845 function plout = getDefaultPlist()
846 persistent pl;
847 if exist('pl', 'var')==0 || isempty(pl)
848 pl = buildplist();
849 end
850 plout = pl;
851 end
852
853 function pl = buildplist()
854 pl = plist({'conf', 'Required percentage confidence level. [0-100]'}, paramValue.DOUBLE_VALUE(95));
855 end
856 % END
857
858 % PARAMETERS:
859 % 'conf' - Required percentage confidence level. It is a
860 % number between 0 and 100 [Default: 95,
861 % corresponding to 95% confidence].