43
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % WELCHSCALE scales the output of welch to be in the required units
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 % $Id: welchscale.m,v 1.1 2011/12/01 09:41:22 hewitson Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 function [yy, dyy, info] = welchscale(xx, dxx, win, fs, norm, inunits)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 nfft = length(win);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 S1 = sum(win);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 S2 = sum(win.^2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 enbw = fs * S2 / (S1*S1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 if isempty(norm)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 norm = 'None';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 switch lower(norm)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 case 'asd'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 yy = sqrt(xx);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 if isempty(dxx)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 dyy = dxx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 dyy = 1./2./sqrt(xx) .* dxx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 info.units = inunits ./ unit('Hz^0.5');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 case 'psd'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 yy = xx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 dyy = dxx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 info.units = inunits.^2/unit('Hz');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 case 'as'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 yy = sqrt(xx * enbw);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 if isempty(dxx)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 dyy = dxx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 dyy = 1./2./sqrt(xx) .* dxx * enbw;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 info.units = inunits;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 case 'ps'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 yy = xx * enbw;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 dyy = dxx * enbw;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 info.units = inunits.^2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 case 'none'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 yy = xx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 dyy = dxx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 info.units = inunits;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 error('Unknown normalisation');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 info.nfft = nfft;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 info.enbw = enbw;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 info.norm = norm;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54
|