0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % MLPSD_M m-file only version 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_M m-file only version of the LPSD algorithm
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 % CALL: [S,Sxx,ENBW] = mlpsd_m(x,f,r,m,L,fs,win,order,olap)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 % VERSION: $Id: mlpsd_m.m,v 1.9 2011/05/22 22:22:50 mauro Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 % HISTORY: 19-05-2007 M Hewitson
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 % Creation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 %
|
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 function varargout = mlpsd_m(varargin)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 import utils.const.*
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 utils.helper.msg(msg.PROC1, 'using MATLAB to compute core DFT');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 % Get inputs
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 x = varargin{1};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 f = varargin{2};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 r = varargin{3};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 m = varargin{4};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 L = varargin{5};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 fs = varargin{6};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 win = varargin{7};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 order = varargin{8};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 olap = varargin{9};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 twopi = 2.0*pi;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 nx = length(x);
|
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
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 disp_each = round(nf/100)*10;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 for jj = 1:nf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 if mod(jj, disp_each) == 0 || jj == 1 || jj == nf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 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
|
46 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 % compute DFT exponent and window
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 l = L(jj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 switch lower(win.type)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 case 'kaiser'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 win = specwin(win.type, l, win.psll);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 win = specwin(win.type, l);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 p = 1i * twopi * m(jj)/l.*[0:l-1];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 C = win.win .* exp(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 Xolap = (1-olap);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 % do segments
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 A = 0.0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 % Compute the number of averages we want here
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 segLen = l;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 nData = length(x);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 ovfact = 1 / (1 - olap);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 davg = (((nData - segLen)) * ovfact) / segLen + 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 navg = round(davg);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 % Compute steps between segments
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 if navg == 1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 shift = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 shift = (nData - segLen) / (navg - 1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 if shift < 1 || isnan(shift)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 shift = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 % disp(sprintf('Seglen: %d\t | Shift: %f\t | navs: %d', segLen, shift, navg))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 start = 1.0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 for ii = 1:navg
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 % compute start index
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 istart = round (start);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 start = start + shift;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 % get segment
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 xs = [x(istart:istart+l-1)].';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 % detrend segment
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 switch order
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 case -1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 % do nothing
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 case 0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 xs = xs - mean(xs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 case 1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 xs = detrend(xs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 xs = polydetrend(xs, order);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 % make DFT
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 a = C*xs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 A = A + a*conj(a);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 if mod(jj, disp_each) == 0 || jj == 1 || jj == nf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 utils.helper.msg(msg.PROC2, 'averaged %d segments', navg);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 A2ns = 2.0*A/navg;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 S1 = win.ws;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 S12 = S1*S1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 S2 = win.ws2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 ENBW(jj) = fs*S2/S12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118 Sxx(jj) = A2ns/fs/S2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 S(jj) = A2ns/S12;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121 end % for j=1:nf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123 varargout{1} = S;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 varargout{2} = Sxx;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 varargout{3} = ENBW;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128
|