0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % MLPSD_MEX calls the ltpda_dft.mex to compute the DFT part of the LPSD algorithm
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 % DESCRIPTION: MLPSD_MEX calls the ltpda_dft.mex to compute the DFT part of the
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % LPSD algorithm
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 % CALL: [S,Sxx,ENBW] = mlpsd_mex(x,f,r,m,L,fs,win,order,olap)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 % VERSION: $Id: mlpsd_mex.m,v 1.13 2011/05/22 22:22:50 mauro Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 % HISTORY: 19-05-2007 M Hewitson
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % Creation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 function varargout = mlpsd_mex(varargin)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 import utils.const.*
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 utils.helper.msg(msg.PROC1, 'using ltpda_dft.mex to compute core DFT');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 % Get inputs
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 x = varargin{1};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 f = varargin{2};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 r = varargin{3};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 m = varargin{4};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 L = varargin{5};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 fs = varargin{6};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 win = varargin{7};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 order = varargin{8};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 olap = varargin{9};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 Lmin = varargin{10};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 twopi = 2.0*pi;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 nf = length(f);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 Sxx = zeros(nf,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 S = zeros(nf,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 ENBW = zeros(nf,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 devxx = zeros(nf,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 dev = zeros(nf,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 disp_each = round(nf/100)*10;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 winType = win.type;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 winPsll = win.psll;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 minReached = 0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 for jj = 1:nf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 if mod(jj, disp_each) == 0 || jj == 1 || jj == nf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 utils.helper.msg(msg.PROC1, 'computing frequency %04d of %d: %f Hz', jj, nf, f(jj));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 % compute DFT exponent and window
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 l = L(jj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 % Check if we need to update the window values
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 % - once we reach Lmin, the window values don't change.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 if ~minReached
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 switch lower(winType)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 case 'kaiser'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 win = specwin(winType, l, winPsll);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 win = specwin(winType, l);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 if l == Lmin
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 minReached = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 p = 1i * twopi * m(jj)/l.*(0:l-1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 C = win.win .* exp(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 % Core DFT part in C-mex file
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 [A, B,nsegs] = ltpda_dft(x, l, C, olap, order);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 if mod(jj, disp_each) == 0 || jj == 1 || jj == nf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 utils.helper.msg(msg.PROC2, 'averaged %d segments', nsegs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 A2ns = 2.0*A;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 B2ns = 4.0*B/nsegs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 S1 = win.ws;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 S12 = S1*S1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 S2 = win.ws2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 ENBW(jj) = fs*S2/S12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 % scale asd/psd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 Sxx(jj) = A2ns/fs/S2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 S(jj) = A2ns/S12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 % scale sqrt(variance)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 devxx(jj) = sqrt(B2ns/fs^2/S2^2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 dev(jj) = sqrt(B2ns/S12^2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 varargout{1} = S;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 varargout{2} = Sxx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 varargout{3} = dev;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 varargout{4} = devxx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 varargout{5} = ENBW;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103
|