Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/spsd.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 % SPSD implements the smoothed (binned) PSD algorithm for analysis objects.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: SPSD implements the smoothed PSD algorithm for analysis objects.
5 %
6 % CALL: bs = spsd(a1,a2,a3,...,pl)
7 % bs = spsd(as,pl)
8 % bs = as.spsd(pl)
9 %
10 % INPUTS: aN - input analysis objects
11 % as - input analysis objects array
12 % pl - input parameter list
13 %
14 % OUTPUTS: bs - array of analysis objects, one for each input
15 %
16 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'spsd')">Parameters Description</a>
17 %
18 % VERSION: $Id: spsd.m,v 1.21 2011/07/11 10:43:35 adrien Exp $
19 %
20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21
22 function varargout = spsd(varargin)
23
24 import utils.const.*
25
26 % Check if this is a call for parameters
27 if utils.helper.isinfocall(varargin{:})
28 varargout{1} = getInfo(varargin{3});
29 return
30 end
31
32 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
33
34 % Collect input variable names
35 in_names = cell(size(varargin));
36 for ii = 1:nargin,in_names{ii} = inputname(ii);end
37
38 % Collect all AOs and plists
39 [as, ao_invars, rest] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
40 [pl, pl_invars, rest] = utils.helper.collect_objects(rest(:), 'plist', in_names);
41
42 % Decide on a deep copy or a modify
43 bs = copy(as, nargout);
44
45 % Combine plists
46 pl = combine(pl, plist(rest(:)), getDefaultPlist);
47
48 inhists = [];
49
50 %% Go through each input AO
51 for jj = 1 : numel(bs)
52 % gather the input history objects
53 inhists = [inhists bs(jj).hist]; %#ok<AGROW>
54
55 % check this is a time-series object
56 if ~isa(bs(jj).data, 'tsdata')
57 warning('!!! spsd requires tsdata (time-series) inputs. Skipping AO %s', ao_invars{jj}); %#ok<WNTAG>
58 else
59
60 % Check the time range.
61 time_range = find(pl, 'times');
62 if ~isempty(time_range)
63 bs(jj) = split(bs(jj), plist('method', 'times', 'times', time_range));
64 end
65 % Check the length of the object
66 if bs(jj).len <= 0
67 error('### The object is empty! Please revise your settings ...');
68 end
69
70 % pl = utils.helper.process_spectral_options(pl, 'log');
71 pl = pl.combine(getDefaultPlist());
72
73 % getting data
74 y = bs(jj).y;
75
76 % Window function
77 Win = find(pl, 'Win');
78 nfft = length(y);
79 Win = ao( combine(plist('win', Win , 'length', nfft), pl) );
80
81 % detrend
82 order = find(pl,'order');
83 if ~(order < 0)
84 y = ltpda_polyreg(y, order).';
85 else
86 y = reshape(y, 1, nfft);
87 end
88
89 % computing PSD
90 window = Win.data.y;
91 window = window/norm(window)*sqrt(nfft);
92 yASD = real(fft(y.*window, nfft)).^2 + imag(fft(y.*window, nfft)).^2;
93 pow = [yASD(1) yASD(2:floor(nfft/2))*2];
94 pow = pow / ( bs(jj).data.fs * nfft);
95 Freqs = linspace(0, bs(jj).data.fs/2, nfft/2);
96
97 % smoothing PSD
98 if ~isempty(find(pl,'frequencies'))
99 error('the option "frequencies" is deprecated, frequencies are "removed" by default')
100 end
101 [Freqs, pow, nFreqs, nDofs] = ltpda_spsd(Freqs, pow, find(pl,'linCoef'), find(pl,'logCoef') );
102 % create new output fsdata
103 scale = find(pl, 'Scale');
104 switch lower(scale)
105 case 'asd'
106 fsd = fsdata(Freqs, sqrt(pow), bs(jj).data.fs);
107 fsd.setYunits(bs(jj).data.yunits / unit('Hz^0.5'));
108 % stdDev = 0.5 * sqrt( pow ./ nDofs ); % linear approximation of the sqrt of a distribution
109 % approximation knowing the STD of the PSD
110 % STD assuming amplitude samples are independent, Chi^1_2 distibuted
111 % (with both variables of powe expectancy pow/2), and of different
112 % magnitude
113 stdDev = 2 * sqrt(pow./nDofs) .* ( nDofs - 2*exp( 2*(gammaln((nDofs+1)/2)-gammaln(nDofs/2)) ) ); % std of the chi_2N^1
114 case 'psd'
115 fsd = fsdata(Freqs, pow, bs(jj).data.fs);
116 fsd.setYunits(bs(jj).data.yunits.^2/unit('Hz'));
117 % STD assuming power samples are independent, Chi^2_2 distibuted
118 % (with both variables of expectancy pow/2), and of different
119 % magnitude
120 stdDev = sqrt(2) * (pow./nDofs) .* sqrt(2*nDofs); % std of the chi_2N^2
121 otherwise
122 error(['### Unknown scaling:' scale]);
123 end
124
125 fsd.setXunits('Hz');
126 fsd.setDx(nFreqs*Freqs(2)/2);
127 fsd.setEnbw(1);% WARNING HERE!!!
128 fsd.setT0(bs(jj).data.t0);
129 % make output analysis object
130 bs(jj).data = fsd;
131 % set name
132 bs(jj).name = ['SPSD(', ao_invars{jj}, ') ' upper(scale)];
133 % Add standard deviation
134 bs(jj).data.dy = stdDev;
135 % Add history
136 bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), inhists(jj));
137
138 end % End tsdata if/else
139 end % End AO loop
140
141 %% Set output
142 if nargout == numel(bs)
143 % List of outputs
144 for ii = 1:numel(bs)
145 varargout{ii} = bs(ii);
146 end
147 else
148 % Single output
149 varargout{1} = bs;
150 end
151
152 end
153
154
155 %--------------------------------------------------------------------------
156 % Get Info Object
157 %--------------------------------------------------------------------------
158 function ii = getInfo(varargin)
159 if nargin == 1 && strcmpi(varargin{1}, 'None')
160 sets = {};
161 pl = [];
162 else
163 sets = {'Default'};
164 pl = getDefaultPlist;
165 end
166 % Build info object
167 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: spsd.m,v 1.21 2011/07/11 10:43:35 adrien Exp $', sets, pl);
168 end
169
170 %--------------------------------------------------------------------------
171 % Get Default Plist
172 %--------------------------------------------------------------------------
173 function pl = getDefaultPlist()
174
175 % Plist for Welch-based, log-scale spaced spectral estimators.
176 pl = plist;
177
178 % Win
179 p = param({'Win',['the window to be applied to the data to remove the ', ...
180 'discontinuities at edges of segments. [default: taken from user prefs] <br>', ...
181 'Only the design parameters of the window object are used. Enter either: <ul>', ...
182 '<li> a specwin window object OR</li>', ...
183 '<li> a string value containing the window name</li></ul>', ...
184 'e.g., <tt>plist(''Win'', ''Kaiser'', ''psll'', 200)</tt>']}, paramValue.WINDOW);
185 pl.append(p);
186
187 % Psll
188 p = param({'Psll',['the peak sidelobe level for Kaiser windows.<br>', ...
189 'Note: it is ignored for all other windows']}, paramValue.DOUBLE_VALUE(200));
190 pl.append(p);
191
192 % Psll
193 p = param({'levelOrder','the contracting order for levelledHanning window'}, paramValue.DOUBLE_VALUE(2));
194 pl.append(p);
195
196 % Order
197 p = param({'Order',['order of segment detrending:<ul>', ...
198 '<li>-1 - no detrending</li>', ...
199 '<li>0 - subtract mean</li>', ...
200 '<li>1 - subtract linear fit</li>', ...
201 '<li>N - subtract fit of polynomial, order N</li></ul>']}, paramValue.DETREND_ORDER);
202 p.val.setValIndex(2);
203 pl.append(p);
204
205 % Times
206 p = param({'Times','time range. If not empty, sets the restricted interval to analyze'}, paramValue.DOUBLE_VALUE([]));
207 pl.append(p);
208
209 % Scale
210 p = param({'Scale',['scaling of output. Choose from:<ul>', ...
211 '<li>PSD - Power Spectral Density</li>', ...
212 '<li>ASD - Amplitude (linear) Spectral Density</li>'...
213 ]}, {1, {'PSD', 'ASD', 'PS', 'AS'}, paramValue.SINGLE});
214 pl.append(p);
215
216 p = param( {'lincoef', 'Linear scale smoothing coefficent (freq. bins)'}, 1);
217 pl.append(p);
218
219 p = param( {'logcoef', ['Logarithmic scale smoothing coefficent<br>', 'Best compromise for both axes is 2/3']}, 2/3);
220 pl.append(p);
221 end