Ticket #2739: ComparePSD.mo

File ComparePSD.mo, 15.8 KB (added by Martin Sjölund, 11 years ago)
Line 
1package ModelicaServices
2 extends Modelica.Icons.Package;
3
4 package Machine
5 extends Modelica.Icons.Package;
6 final constant Real eps = 1.e-15;
7 final constant Real small = 1.e-60;
8 final constant Real inf = 1.e+60;
9 final constant Integer Integer_inf = OpenModelica.Internal.Architecture.integerMax();
10 end Machine;
11end ModelicaServices;
12
13package Modelica
14 extends Modelica.Icons.Package;
15
16 package Blocks
17 extends Modelica.Icons.Package;
18
19 package Interfaces
20 extends Modelica.Icons.InterfacesPackage;
21 connector RealOutput = output Real;
22
23 partial block SO
24 extends Modelica.Blocks.Icons.Block;
25 RealOutput y;
26 end SO;
27 end Interfaces;
28
29 package Icons
30 extends Modelica.Icons.IconsPackage;
31
32 partial block Block end Block;
33 end Icons;
34 end Blocks;
35
36 package Math
37 extends Modelica.Icons.Package;
38
39 package Icons
40 extends Modelica.Icons.IconsPackage;
41
42 partial function AxisCenter end AxisCenter;
43 end Icons;
44
45 function asin
46 extends Modelica.Math.Icons.AxisCenter;
47 input Real u;
48 output .Modelica.SIunits.Angle y;
49 external "builtin" y = asin(u);
50 end asin;
51
52 function exp
53 extends Modelica.Math.Icons.AxisCenter;
54 input Real u;
55 output Real y;
56 external "builtin" y = exp(u);
57 end exp;
58 end Math;
59
60 package Constants
61 extends Modelica.Icons.Package;
62 final constant Real pi = 2 * Math.asin(1.0);
63 final constant .Modelica.SIunits.Velocity c = 299792458;
64 final constant Real mue_0(final unit = "N/A2") = 4 * pi * 1.e-7;
65 end Constants;
66
67 package Icons
68 extends Icons.Package;
69
70 partial package ExamplesPackage
71 extends Modelica.Icons.Package;
72 end ExamplesPackage;
73
74 partial model Example end Example;
75
76 partial package Package end Package;
77
78 partial package InterfacesPackage
79 extends Modelica.Icons.Package;
80 end InterfacesPackage;
81
82 partial package UtilitiesPackage
83 extends Modelica.Icons.Package;
84 end UtilitiesPackage;
85
86 partial package IconsPackage
87 extends Modelica.Icons.Package;
88 end IconsPackage;
89 end Icons;
90
91 package SIunits
92 extends Modelica.Icons.Package;
93
94 package Conversions
95 extends Modelica.Icons.Package;
96
97 package NonSIunits
98 extends Modelica.Icons.Package;
99 type Temperature_degC = Real(final quantity = "ThermodynamicTemperature", final unit = "degC");
100 end NonSIunits;
101 end Conversions;
102
103 type Angle = Real(final quantity = "Angle", final unit = "rad", displayUnit = "deg");
104 type Time = Real(final quantity = "Time", final unit = "s");
105 type Velocity = Real(final quantity = "Velocity", final unit = "m/s");
106 type Acceleration = Real(final quantity = "Acceleration", final unit = "m/s2");
107 type Period = Real(final quantity = "Time", final unit = "s");
108 type Frequency = Real(final quantity = "Frequency", final unit = "Hz");
109 type FaradayConstant = Real(final quantity = "FaradayConstant", final unit = "C/mol");
110 end SIunits;
111end Modelica;
112
113package Noise
114 extends Modelica.Icons.Package;
115
116 model GlobalSeed
117 parameter Integer userSeed = 1;
118 final parameter Integer seed = userSeed;
119 end GlobalSeed;
120
121 block PRNG
122 extends Modelica.Blocks.Interfaces.SO;
123 outer GlobalSeed globalSeed;
124 parameter Boolean useSampleBasedMethods = false;
125 replaceable function SampleBasedRNG = Noise.RNG.SampleBased.RNG_LCG constrainedby Noise.Utilities.Interfaces.SampleBasedRNG;
126 replaceable function SampleFreeRNG = Noise.RNG.SampleFree.RNG_DIRCS constrainedby Noise.Utilities.Interfaces.SampleFreeRNG;
127 protected
128 function SampleBasedRNG0 = SampleBasedRNG;
129 function SampleFreeRNG0 = SampleFreeRNG;
130 public
131 replaceable function PDF = Noise.PDF.PDF_Uniform constrainedby Noise.Utilities.Interfaces.PDF;
132 protected
133 function SampleBasedPDF0 = PDF(redeclare function RNG = SampleBasedRNG0);
134 function SampleFreePDF0 = PDF(redeclare function RNG = SampleFreeRNG0);
135 public
136 parameter Boolean infiniteFreq = false;
137 protected
138 parameter Modelica.SIunits.Frequency freq = 0.5 * 1 / samplePeriod;
139 public
140 replaceable function PSD = Noise.PSD.PSD_WhiteNoise constrainedby Noise.Utilities.Interfaces.PSD;
141 protected
142 function SampleBasedPSD0 = PSD(redeclare function PDF = SampleBasedPDF0);
143 function SampleFreePSD0 = PSD(redeclare function PDF = SampleFreePDF0);
144 function InfiniteFreqPSD0 = Noise.PSD.PSD_WhiteNoise(redeclare function PDF = SampleFreePDF0);
145 public
146 parameter Modelica.SIunits.Time startTime = 0;
147 parameter Modelica.SIunits.Time samplePeriod = 0.01;
148 parameter Boolean enable = true;
149 parameter Real y_off = 0;
150 replaceable function Seed = Noise.Seed.Seed_MRG(real_seed = 0.0) constrainedby Noise.Utilities.Interfaces.Seed;
151 protected
152 parameter Integer state_size = 33;
153 Integer[state_size] state;
154 Real t_last;
155 public
156 parameter Integer localSeed = 123456789;
157 parameter Boolean useGlobalSeed = true;
158 final parameter Integer seed = if useGlobalSeed then Utilities.Auxiliary.combineSeedLCG(localSeed, globalSeed.seed) else localSeed;
159 final parameter Real DT = 1 / (2 * freq);
160 output Real y_hold;
161 protected
162 discrete Real dummy1;
163 discrete Real dummy2;
164 initial equation
165 if useSampleBasedMethods then
166 pre(state) = Seed(local_seed = localSeed, global_seed = if useGlobalSeed then globalSeed.seed else 0, n = state_size, real_seed = 0.0);
167 pre(t_last) = floor(time / DT) * DT;
168 end if;
169 equation
170 if not enable then
171 y = y_off;
172 y_hold = y_off;
173 t_last = 0;
174 dummy1 = 0;
175 dummy2 = 0;
176 state = zeros(state_size);
177 else
178 if useSampleBasedMethods then
179 when sample(0, DT) then
180 t_last = time;
181 (dummy1, dummy2, state) = SampleBasedPSD0(t = time, dt = DT, t_last = pre(t_last), states_in = pre(state));
182 end when;
183 (y_hold, y) = SampleBasedPSD0(t = time, dt = DT, t_last = t_last, states_in = state);
184 else
185 when initial() then
186 dummy1 = 0;
187 dummy2 = 0;
188 end when;
189 state = Seed(local_seed = localSeed, global_seed = if useGlobalSeed then globalSeed.seed else 0, n = state_size, real_seed = 0.0);
190 t_last = noEvent(2 * abs(time) + 1);
191 if infiniteFreq then
192 (y_hold, y) = InfiniteFreqPSD0(t = time, dt = 0, t_last = t_last, states_in = state);
193 else
194 (y_hold, y) = SampleFreePSD0(t = time, dt = DT, t_last = t_last, states_in = state);
195 end if;
196 end if;
197 end if;
198 end PRNG;
199
200 package RNG
201 extends Modelica.Icons.Package;
202
203 package SampleBased
204 extends Modelica.Icons.Package;
205
206 function RNG_MRG
207 extends Noise.Utilities.Interfaces.SampleBasedRNG;
208 input Integer[:] a = {1071064, 0, 0, 0, 0, 0, 2113664};
209 input Integer c = 0;
210 input Integer m = 2147483629;
211 algorithm
212 assert(size(states_in, 1) >= size(a, 1), "State must have at least as many elements as a!");
213 states_out := states_in;
214 states_out[1] := 0;
215 for i in 1:size(a, 1) loop
216 states_out[1] := states_out[1] + a[i] * states_in[i];
217 end for;
218 states_out[1] := integer(mod(states_out[1] + c, m));
219 for i in 1:size(a, 1) - 1 loop
220 states_out[i + 1] := states_in[i];
221 end for;
222 rand := abs(states_out[1] / (m - 1));
223 end RNG_MRG;
224
225 function RNG_LCG
226 extends Noise.Utilities.Interfaces.SampleBasedRNG;
227 input Integer a = 69069;
228 input Integer c = 1;
229 input Integer m = 2147483647;
230 algorithm
231 (rand, states_out) := RNG_MRG(instance, states_in, a = {a}, c = c, m = m);
232 end RNG_LCG;
233 end SampleBased;
234
235 package SampleFree
236 extends Modelica.Icons.Package;
237
238 function RNG_DIRCS
239 extends Noise.Utilities.Interfaces.SampleFreeRNG;
240 replaceable function Seed = Noise.Seed.Seed_Real constrainedby Noise.Utilities.Interfaces.Seed;
241 replaceable function RNG = Noise.RNG.SampleBased.RNG_MRG(a = {134775813, 134775813}, c = 1) constrainedby Noise.Utilities.Interfaces.RNG;
242 input Integer k = 1;
243 protected
244 Integer[2] states_internal;
245 algorithm
246 states_internal := Seed(real_seed = instance, local_seed = states_in[1], global_seed = 0, n = 2);
247 for i in 1:k loop
248 (rand, states_internal) := RNG(instance = instance, states_in = states_internal);
249 end for;
250 states_out := states_in;
251 end RNG_DIRCS;
252 end SampleFree;
253 end RNG;
254
255 package PDF
256 extends Noise.Utilities.Icons.PDFPackage;
257
258 function PDF_Uniform
259 extends Noise.Utilities.Interfaces.PDF;
260 input Real[2] interval = {0, 1};
261 algorithm
262 (rand, states_out) := RNG(instance = instance, states_in = states_in);
263 rand := rand * (interval[2] - interval[1]) + interval[1];
264 end PDF_Uniform;
265 end PDF;
266
267 package PSD
268 extends Noise.Utilities.Icons.PSDPackage;
269
270 function PSD_WhiteNoise
271 extends Noise.Utilities.Interfaces.PSD;
272 algorithm
273 if dt > 0 then
274 (rand, states_out) := PDF(instance = floor(t / dt) * dt, states_in = states_in);
275 else
276 (rand, states_out) := PDF(instance = t, states_in = states_in);
277 end if;
278 rand_hold := rand;
279 end PSD_WhiteNoise;
280
281 function PSD_IdealLowPass
282 extends PSD_Interpolation(redeclare function Kernel = Kernels.IdealLowPass);
283 end PSD_IdealLowPass;
284
285 function PSD_LinearInterpolation
286 extends PSD_Interpolation(redeclare function Kernel = Kernels.Linear, n = 1);
287 end PSD_LinearInterpolation;
288
289 function PSD_Interpolation
290 extends Noise.Utilities.Interfaces.PSD;
291 replaceable function Kernel = Noise.PSD.Kernels.IdealLowPass constrainedby Utilities.Interfaces.Kernel;
292 input Integer n = 5;
293 input Integer max_n = n;
294 protected
295 Real raw;
296 Real coefficient;
297 Real scaling;
298 Integer[size(states_in, 1)] states_temp;
299 algorithm
300 rand := 0;
301 scaling := 0;
302 states_temp := states_in;
303 for i in (-max_n):(-n) loop
304 (raw, states_temp) := PDF(instance = (floor(t / dt) + i) * dt, states_in = states_temp);
305 end for;
306 for i in (-n) + 1:n loop
307 (raw, states_temp) := PDF(states_in = states_temp, instance = floor(t / dt + i) * dt);
308 coefficient := if t_last <= t then Kernel(t = t - (t_last + i * dt), dt = dt) else Kernel(t = t - floor(t / dt + i) * dt, dt = dt);
309 rand := rand + raw * coefficient;
310 scaling := scaling + coefficient;
311 if i == 0 then
312 rand_hold := raw;
313 else
314 end if;
315 end for;
316 rand := rand / scaling;
317 (raw, states_out) := PDF(states_in = states_in, instance = floor(t / dt) * dt);
318 end PSD_Interpolation;
319
320 package Kernels
321 extends Modelica.Icons.Package;
322
323 function IdealLowPass
324 extends Noise.Utilities.Interfaces.Kernel;
325 input Modelica.SIunits.Frequency B = 1 / 2 / dt;
326 algorithm
327 h := 2 * B * .Noise.Utilities.Math.sinc(2 * .Modelica.Constants.pi * B * t);
328 end IdealLowPass;
329
330 function Linear
331 extends Noise.Utilities.Interfaces.Kernel;
332 algorithm
333 h := if t < (-dt) then 0 else if t < 0 then 1 + t / dt else if t < dt then 1 - t / dt else 0;
334 end Linear;
335 end Kernels;
336 end PSD;
337
338 package Seed
339 extends Noise.Utilities.Icons.SeedPackage;
340
341 function Seed_MRG
342 extends Utilities.Interfaces.Seed;
343 input Integer[:] a = fill(134775813, n);
344 input Integer c = 1;
345 input Integer m = 2147483629;
346 input Integer k = n;
347 protected
348 Real dummy;
349 Integer[max(n, 2)] internal_states;
350 algorithm
351 assert(n > 0, "You are seeding a state vector of size 0!");
352 internal_states := cat(1, {local_seed, global_seed}, fill(0, max(n, 2) - 2));
353 for i in 1:k loop
354 (dummy, internal_states) := RNG.SampleBased.RNG_MRG(instance = real_seed, states_in = internal_states, a = a, c = c, m = m);
355 end for;
356 for i in 1:n loop
357 states[i] := internal_states[i];
358 end for;
359 end Seed_MRG;
360
361 function Seed_Real
362 extends Utilities.Interfaces.Seed;
363 algorithm
364 states := Noise.Utilities.Auxiliary.SeedReal(local_seed = local_seed, global_seed = global_seed, real_seed = real_seed, n = n);
365 end Seed_Real;
366 end Seed;
367
368 package Utilities
369 extends Modelica.Icons.Package;
370 extends Modelica.Icons.UtilitiesPackage;
371
372 package Icons
373 extends Modelica.Icons.IconsPackage;
374
375 partial function PDF end PDF;
376
377 partial package PDFPackage
378 extends Modelica.Icons.Package;
379 end PDFPackage;
380
381 partial function PSD end PSD;
382
383 partial package PSDPackage
384 extends Modelica.Icons.Package;
385 end PSDPackage;
386
387 partial function Seed end Seed;
388
389 partial package SeedPackage
390 extends Modelica.Icons.Package;
391 end SeedPackage;
392 end Icons;
393
394 package Interfaces
395 extends Modelica.Icons.InterfacesPackage;
396
397 partial function InputOutput
398 input Modelica.SIunits.Time instance;
399 input Integer[:] states_in;
400 output Real rand;
401 output Integer[size(states_in, 1)] states_out;
402 end InputOutput;
403
404 partial function RNG
405 extends Interfaces.InputOutput;
406 end RNG;
407
408 partial function SampleBasedRNG
409 extends Interfaces.RNG;
410 end SampleBasedRNG;
411
412 partial function SampleFreeRNG
413 extends Interfaces.RNG;
414 end SampleFreeRNG;
415
416 partial function PDF
417 extends Icons.PDF;
418 extends Interfaces.InputOutput;
419 replaceable function RNG = Noise.RNG.SampleBased.RNG_LCG constrainedby Interfaces.RNG;
420 end PDF;
421
422 partial function PSD
423 extends Icons.PSD;
424 output Real rand_hold;
425 extends Interfaces.InputOutput(instance = t);
426 input Modelica.SIunits.Time t;
427 input Modelica.SIunits.Period dt;
428 input Modelica.SIunits.Time t_last;
429 replaceable function PDF = Noise.PDF.PDF_Uniform constrainedby Interfaces.PDF;
430 end PSD;
431
432 partial function Kernel
433 input Real t;
434 input Real dt;
435 output Real h;
436 end Kernel;
437
438 partial function Seed
439 extends Icons.Seed;
440 input Integer local_seed = 12345;
441 input Integer global_seed = 67890;
442 input Real real_seed = 1.234;
443 input Integer n = 33;
444 output Integer[n] states;
445 end Seed;
446
447 partial function combineSeed
448 input Integer seed1;
449 input Integer seed2;
450 output Integer newSeed;
451 end combineSeed;
452 end Interfaces;
453
454 package Auxiliary
455 extends Modelica.Icons.Package;
456
457 function SeedReal
458 input Integer local_seed;
459 input Integer global_seed;
460 input Real real_seed;
461 input Integer n;
462 output Integer[n] states;
463 external "C" NOISE_SeedReal(local_seed, global_seed, real_seed, n, states);
464 end SeedReal;
465
466 function combineSeedLCG
467 extends Interfaces.combineSeed;
468 external "C" newSeed = NOISE_combineSeedLCG(seed1, seed2);
469 end combineSeedLCG;
470 end Auxiliary;
471
472 package Math
473 extends Modelica.Icons.Package;
474
475 function sinc
476 input Real x;
477 output Real y;
478 algorithm
479 y := if abs(x) > 0.5e-4 then sin(x) / x else 1 - x ^ 2 / 6 + x ^ 4 / 120;
480 end sinc;
481 end Math;
482 end Utilities;
483end Noise;
484
485model ComparePSD
486 extends Modelica.Icons.Example;
487 .Noise.PRNG WhiteNoise(redeclare function PSD = .Noise.PSD.PSD_WhiteNoise, useSampleBasedMethods = false, redeclare function PDF = .Noise.PDF.PDF_Uniform(interval = {-1, 1}));
488 .Noise.PRNG IdealLowPass(redeclare function PSD = .Noise.PSD.PSD_IdealLowPass(n = 10), useSampleBasedMethods = false, redeclare function PDF = .Noise.PDF.PDF_Uniform(interval = {-1, 1}));
489 .Noise.PRNG Linear(redeclare function PSD = .Noise.PSD.PSD_LinearInterpolation(n = 5), useSampleBasedMethods = false, redeclare function PDF = .Noise.PDF.PDF_Uniform(interval = {-1, 1}));
490 inner .Noise.GlobalSeed globalSeed;
491end ComparePSD;