Mercurial > hg > ltpda

comparison m-toolbox/classes/@ssm/noiseSpectrum.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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %
3 % DESCRIPTION: noiseSpectrum makes spectra of ssm outputs from the given noises spectra.
4 %
5 % CALL: [ao_out] = noiseSpectrum(sys, plist_inputs)
6 %
7 % INPUTS:
8 % - sys, (array of) ssm object
9 % - plist_inputs contains names of noise inputs, spectra of
10 % specified noise inputs, returned
11 % outputs and frequency's vector of returned spectrum
12 %
13 % OUTPUTS:
14 % _ ao_out contains spectra of specified outputs and states
15 %
16 % <a href="matlab:utils.helper.displayMethodInfo('ssm', 'noiseSpectrum')">Parameters Description</a>
17 %
18 % VERSION: $Id: noiseSpectrum.m,v 1.8 2011/04/08 08:56:23 hewitson Exp $
19 %
20 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21
22
23
24 function varargout = noiseSpectrum(varargin)
25
26 %% starting initial checks
27 error('### this function is deprecated, use ssm/PSD or ssm/CPSD instead')
28
29 % Check if this is a call for parameters
30 if utils.helper.isinfocall(varargin{:})
31 varargout{1} = getInfo(varargin{3});
32 return
33 end
34
35 utils.helper.msg(utils.const.msg.MNAME, ['running ', mfilename]);
36
37 % Collect input variable names
38 in_names = cell(size(varargin));
39 for ii = 1:nargin,in_names{ii} = inputname(ii);end
40
41 % Collect all SSMs and plists
42 [sys, ssm_invars, rest] = utils.helper.collect_objects(varargin(:), 'ssm', in_names);
43 [pl, invars2, rest] = utils.helper.collect_objects(rest(:), 'plist');
44 pl = combine(pl, plist(rest{:}) );
45 plist_inputs = combine(pl, getDefaultPlist());
46
47
48 %% retrieve system infos
49 if numel(sys)~=1
50 error('noisespectrum needs exactly one ssm as an input')
51 end
52 if ~sys.isnumerical
53 error(['error because system ',sys.name,' is not numerical']);
54 end
55
56 %% collecting data
57 noise_varnames = find(plist_inputs, 'noise variable names');
58 if ischar(noise_varnames)
59 noise_varnames = {noise_varnames};
60 end
61 Nnoise = numel(noise_varnames);
62 In_noise_spectrum = find(plist_inputs, 'spectrum');
63 freq_vector = find(plist_inputs, 'f');
64 return_outputs = find(plist_inputs, 'return outputs');
65 if ischar(return_outputs)
66 return_outputs = {return_outputs};
67 end
68 Noutputs=numel(return_outputs);
69
70 %% verify frequency
71 Nspectrum=length(In_noise_spectrum);
72 ur=1;
73 diff_vec=[];
74 vec=[];
75 for u=1: Nspectrum
76 noise=In_noise_spectrum{u};
77 if isa(noise,'ao')==1
78 vec(ur,:)=noise.x;
79 ur=ur+1;
80 end
81 end
82
83 if ~isempty(vec) && ur>=3
84 for ui=1: ur-2
85 diff_vec(ui,:)= vec(ui+1,:)-vec(ui,:);
86 end
87 end
88
89 if ~isempty(diff_vec)
90 if any(diff_vec~=0)
91 freq_status=0;
92 error('### Please specify inputs in the same frequency range .')
93 else
94 freq_status=1;
95 end
96 end
97
98 if isempty(vec) && isempty(freq_vector)
99 f1 = find(plist_inputs,'f1');
100 f2 = find(plist_inputs,'f2');
101 nf = find(plist_inputs,'nf');
102 if isempty(f1) || isempty(f2)|| isempty(nf)
103 error('### Please specify frequency vector a start and stop frequency .');
104 else
105 freq_vector = linspace(f1, f2, nf);
106 end
107
108 end
109
110 if ~isempty(freq_vector)
111 if ur==2 || (ur>=3 && freq_status==1)
112 if length(freq_vector) ~=length(vec(1,:)) || any((freq_vector - vec(1,:))~=0)
113 error('### Please specify inputs in asked frequency vector .');
114 else
115 freq_vector=vec(1,:);
116 end
117 end
118 else
119 if ur==2 || (ur>=3 && freq_status==1)
120 freq_vector=vec(1,:);
121 end
122 end
123
124 %% spectra inputs
125 if Nnoise>=1
126 for u=1: Nnoise
127 noise=In_noise_spectrum{u};
128 if isa(noise,'ao')==1
129 noise_uncor{u}= In_noise_spectrum{u}.y';
130 elseif isa(noise,'double')==1
131 noise_uncor{u}=In_noise_spectrum{u}.*ones(1,length(freq_vector));
132 elseif isa(noise,'miir')==1
133 rep=resp(In_noise_spectrum{u},plist('f',freq_vector));
134 noise_uncor{u}= rep.y';
135 elseif isa(noise,'pzmodel')==1
136 rep=resp(In_noise_spectrum{u},plist('f',freq_vector));
137 noise_uncor{u}=rep.y';
138 else error('spectrum should be an Ao, a double, a miir or a pzmodel')
139 end
140 end
141 else
142 error('Specify at least one input noise spectrum')
143 end
144
145 %% SSM Transfer function
146 plist_OutputTF=plist('inputs', noise_varnames,'outputs', return_outputs, 'f',freq_vector);
147 if ~isempty(ver('control')) % check that control system toolbox is installed
148 rep_TF= sys.bodecst(plist_OutputTF);
149 else
150 rep_TF= sys.bode(plist_OutputTF);
151 end
152
153 [Ninp Nout]=size(rep_TF);
154 if (Ninp~=Nnoise)||(Nout~=Noutputs)
155 error('Asked input /output does not exist. Please verify names of return outputs and noise variable names')
156 end
157
158 %% Spectral density
159 dd=[noise_uncor{1,:}];
160 noise=reshape(dd,length(freq_vector),Nnoise);
161 for u=1:length(freq_vector)
162 mat_TF=reshape(rep_TF.y(u),Noutputs,Nnoise);
163 for ui=1:Noutputs
164 mat_TFN(ui,:)= mat_TF(ui,:).* abs(noise(u,:));
165 end
166 resu(:,:,u)=mat_TFN*conj(transpose(mat_TF));
167 end
168
169 %% saving in aos
170
171 ao_out = ao.initObjectWithSize(Noutputs,Noutputs);
172 for i=1:Noutputs
173
174 for j=1:Noutputs
175 if j==i
176 nameao=[sprintf('PSD(%s)', return_outputs{i})];
177 else
178 nameao=[sprintf('CPD(%s)', [return_outputs{i} ',' return_outputs{j}])];
179 end
180 for u=1:length(freq_vector)
181 ao_out_y(u)= resu(j,i,u);
182 end
183 ao_out(i,j) = ao(fsdata(freq_vector,ao_out_y));
184 ao_out(i,j).setName(nameao);
185 ao_out(i,j).setXunits('Hz');
186 ao_out(i,j).setYunits('arb');
187 ao_out(i,j).setDescription(...
188 ['Cross Power Spectrum of ' , return_outputs{j} , return_outputs{i}]);
189 end
190 if nargout == numel(ao_out)
191 for ii = 1:numel(ao_out)
192 varargout{ii} = ao_out(ii);
193 end
194 else
195 varargout{1} = ao_out;
196 end
197 end
198 varargout = {ao_out };
199 end
200
201
202 %--------------------------------------------------------------------------
203 % Get Info Object
204 %--------------------------------------------------------------------------
205 function ii = getInfo(varargin)
206
207 if nargin == 1 && strcmpi(varargin{1}, 'None')
208 sets = {};
209 pl = [];
210 else
211 sets = {'Default'};
212 pl = getDefaultPlist;
213 end
214 % Build info object
215 ii = minfo(mfilename, 'ssm', 'ltpda', utils.const.categories.op, '$Id: noiseSpectrum.m,v 1.8 2011/04/08 08:56:23 hewitson Exp $', sets, pl);
216
217 end
218
219 %--------------------------------------------------------------------------
220 % Get Default Plist
221 %--------------------------------------------------------------------------
222 function pl = getDefaultPlist()
223 pl = plist();
224
225 p = param({'noise variable names', 'A cell-array of strings specifying the desired input variable names.'}, {} );
226 pl.append(p);
227
228 p = param({'spectrum', 'The spectrum of this noise between input ports for the <i>time-continuous</i> noise model.'}, []);
229 pl.append(p);
230
231 p = param({'f', 'frequency vector of output spectrum .'}, []);
232 pl.append(p);
233
234 p = param({'return outputs', 'A cell-array of output ports to return.'}, {});
235 pl.append(p);
236
237 end
238