Ruprecht Altenburger / Cntrlol_Lib

My Controller Library

GPA.cpp@1:bf62e74fbcf3, 2018-09-28 (annotated)

Committer:
altb
Date:
Fri Sep 28 09:01:42 2018 +0000
Revision:
1:bf62e74fbcf3
Parent:
0:e2a7d7f91e49 
.

Who changed what in which revision?

UserRevisionLine numberNew contents of line
altb 0:e2a7d7f91e49 1 /*
altb 0:e2a7d7f91e49 2 GPA: frequency point wise gain and phase analyser to measure the frequency
altb 0:e2a7d7f91e49 3 respone of a dynamical system
altb 0:e2a7d7f91e49 4
altb 0:e2a7d7f91e49 5 hint: the measurements should only be perfomed in closed loop
altb 0:e2a7d7f91e49 6 assumption: the system is at the desired steady state of interest when
altb 0:e2a7d7f91e49 7 the measurment starts
altb 0:e2a7d7f91e49 8
altb 0:e2a7d7f91e49 9 exc(2) C: controller
altb 0:e2a7d7f91e49 10 | P: plant
altb 0:e2a7d7f91e49 11 v
altb 0:e2a7d7f91e49 12 exc(1) --> o ->| C |--->o------->| P |----------> out
altb 0:e2a7d7f91e49 13 ^ - | |
altb 0:e2a7d7f91e49 14 | --> inp | exc: excitation signal (E)
altb 0:e2a7d7f91e49 15 | | inp: input plant (U)
altb 0:e2a7d7f91e49 16 -------------------------------- out: output plant (Y)
altb 0:e2a7d7f91e49 17
altb 0:e2a7d7f91e49 18 instantiate option 1: (logarithmic equaly spaced frequency points)
altb 0:e2a7d7f91e49 19
altb 0:e2a7d7f91e49 20 GPA(float fMin, float fMax, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1)
altb 0:e2a7d7f91e49 21
altb 0:e2a7d7f91e49 22 fMin: minimal desired frequency that should be measured in Hz
altb 0:e2a7d7f91e49 23 fMax: maximal desired frequency that should be measured in Hz
altb 0:e2a7d7f91e49 24 NfexcDes: number of logarithmic equaly spaced frequency points
altb 0:e2a7d7f91e49 25 NperMin: minimal number of periods that are used for each frequency point
altb 0:e2a7d7f91e49 26 NmeasMin: minimal number of samples that are used for each frequency point
altb 0:e2a7d7f91e49 27 Ts: sampling time
altb 0:e2a7d7f91e49 28 Aexc0: excitation amplitude at fMin
altb 0:e2a7d7f91e49 29 Aexc1: excitation amplitude at fMax
altb 0:e2a7d7f91e49 30
altb 0:e2a7d7f91e49 31 instantiate option 2: (for a second, refined frequency grid measurement)
altb 0:e2a7d7f91e49 32
altb 0:e2a7d7f91e49 33 GPA(float f0, float f1, float *fexcDes, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1)
altb 0:e2a7d7f91e49 34
altb 0:e2a7d7f91e49 35 f0: frequency point for the calculation of Aexc0 in Hz (should be chosen like in the first measurement)
altb 0:e2a7d7f91e49 36 f1: frequency point for the calculation of Aexc1 in Hz (should be chosen like in the first measurement)
altb 0:e2a7d7f91e49 37 *fexcDes: sorted frequency point array in Hz
altb 0:e2a7d7f91e49 38 NfexcDes: length of fexcDes
altb 0:e2a7d7f91e49 39
altb 0:e2a7d7f91e49 40 instantiate option 3: (for an arbitary but sorted frequency grid measurement)
altb 0:e2a7d7f91e49 41
altb 0:e2a7d7f91e49 42 GPA(float *fexcDes, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1)
altb 0:e2a7d7f91e49 43
altb 0:e2a7d7f91e49 44 *fexcDes: sorted frequency point array in Hz
altb 0:e2a7d7f91e49 45 Aexc0: excitation amplitude at fexcDes[0]
altb 0:e2a7d7f91e49 46 Aexc1: excitation amplitude at fexcDes[NfexcDes-1]
altb 0:e2a7d7f91e49 47 NfexcDes: length of fexcDes
altb 0:e2a7d7f91e49 48
altb 0:e2a7d7f91e49 49 hints: the amplitude drops with 1/fexc, if you're using Axc1 = Aexc0/fMax then d/dt exc = const.,
altb 0:e2a7d7f91e49 50 this is recommended if your controller does not have a rolloff. if a desired frequency point
altb 0:e2a7d7f91e49 51 is not measured try to increase Nmeas.
altb 0:e2a7d7f91e49 52
altb 0:e2a7d7f91e49 53 pseudo code for a closed loop measurement with a controller C:
altb 0:e2a7d7f91e49 54
altb 0:e2a7d7f91e49 55 excitation input at (1):
altb 0:e2a7d7f91e49 56
altb 0:e2a7d7f91e49 57 - measuring the plant P and the closed loop tf T = PC/(1 + PC):
altb 0:e2a7d7f91e49 58 desTorque = pi_w(omega_desired - omega + excWobble);
altb 0:e2a7d7f91e49 59 inpWobble = desTorque;
altb 0:e2a7d7f91e49 60 outWobble = omega;
altb 0:e2a7d7f91e49 61 excWobble = Wobble(excWobble, outWobble);
altb 0:e2a7d7f91e49 62
altb 0:e2a7d7f91e49 63 - measuring the controller C and the closed loop tf SC = C/(1 + PC)
altb 0:e2a7d7f91e49 64 desTorque = pi_w(omega_desired - omega + excWobble);
altb 0:e2a7d7f91e49 65 inpWobble = n_soll + excWobble - omega;
altb 0:e2a7d7f91e49 66 outWobble = desTorque;
altb 0:e2a7d7f91e49 67 excWobble = Wobble(inpWobble, outWobble);
altb 0:e2a7d7f91e49 68
altb 0:e2a7d7f91e49 69 excitation input at (2):
altb 0:e2a7d7f91e49 70
altb 0:e2a7d7f91e49 71 - measuring the plant P and the closed loop tf SP = P/(1 + PC):
altb 0:e2a7d7f91e49 72 desTorque = pi_w(omega_desired - omega) + excWobble;
altb 0:e2a7d7f91e49 73 inpWobble = desTorque;
altb 0:e2a7d7f91e49 74 outWobble = omega;
altb 0:e2a7d7f91e49 75 excWobble = Wobble(excWobble, outWobble);
altb 0:e2a7d7f91e49 76
altb 0:e2a7d7f91e49 77 usage:
altb 0:e2a7d7f91e49 78 exc(k+1) = myGPA(inp(k), out(k)) does update the internal states of the
altb 0:e2a7d7f91e49 79 gpa at the timestep k and returns the excitation signal for the timestep
altb 0:e2a7d7f91e49 80 k+1. the results are plotted to a terminal (putty) over a serial
altb 0:e2a7d7f91e49 81 connection and look as follows:
altb 0:e2a7d7f91e49 82 -----------------------------------------------------------------------------------------
altb 0:e2a7d7f91e49 83 fexc[Hz] |Gyu| ang(Gyu) |Gye| ang(Gye) |E| |U| |Y|
altb 0:e2a7d7f91e49 84 -----------------------------------------------------------------------------------------
altb 0:e2a7d7f91e49 85 1.000e+00 2.924e+01 -1.572e+00 1.017e+00 -4.983e-02 5.000e+00 1.739e-01 5.084e+00
altb 0:e2a7d7f91e49 86
altb 0:e2a7d7f91e49 87 in matlab you can use:
altb 0:e2a7d7f91e49 88 dataNucleo = [... insert measurement data here ...];
altb 0:e2a7d7f91e49 89 g = frd(dataNucleo(:,2).*exp(1i*dataNucleo(:,3)), dataNucleo(:,1), Ts, 'Units', 'Hz');
altb 0:e2a7d7f91e49 90 gcl = frd(dataNucleo(:,4).*exp(1i*dataNucleo(:,5)), dataNucleo(:,1), Ts, 'Units', 'Hz');
altb 0:e2a7d7f91e49 91
altb 0:e2a7d7f91e49 92 if you're evaluating more than one measurement which contain equal frequency points try:
altb 0:e2a7d7f91e49 93 dataNucleo = [dataNucleo1; dataNucleo2];
altb 0:e2a7d7f91e49 94 [~, ind] = unique(dataNucleo(:,1),'stable');
altb 0:e2a7d7f91e49 95 dataNucleo = dataNucleo(ind,:);
altb 0:e2a7d7f91e49 96
altb 0:e2a7d7f91e49 97 autor: M.E. Peter
altb 0:e2a7d7f91e49 98 */
altb 0:e2a7d7f91e49 99
altb 0:e2a7d7f91e49 100 #include "GPA.h"
altb 0:e2a7d7f91e49 101 #include "mbed.h"
altb 0:e2a7d7f91e49 102 #include "math.h"
altb 0:e2a7d7f91e49 103 #define pi 3.141592653589793
altb 0:e2a7d7f91e49 104
altb 0:e2a7d7f91e49 105 using namespace std;
altb 0:e2a7d7f91e49 106
altb 0:e2a7d7f91e49 107 GPA::GPA(float fMin, float fMax, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1)
altb 0:e2a7d7f91e49 108 {
altb 0:e2a7d7f91e49 109 this->NfexcDes = NfexcDes;
altb 0:e2a7d7f91e49 110 this->NperMin = NperMin;
altb 0:e2a7d7f91e49 111 this->NmeasMin = NmeasMin;
altb 0:e2a7d7f91e49 112 this->Ts = (double)Ts;
altb 0:e2a7d7f91e49 113
altb 0:e2a7d7f91e49 114 // calculate logarithmic spaced frequency points
altb 0:e2a7d7f91e49 115 fexcDes = (double*)malloc(NfexcDes*sizeof(double));
altb 0:e2a7d7f91e49 116 fexcDesLogspace((double)fMin, (double)fMax, NfexcDes);
altb 0:e2a7d7f91e49 117
altb 0:e2a7d7f91e49 118 // calculate coefficients for decreasing amplitude (1/fexc)
altb 0:e2a7d7f91e49 119 this->aAexcDes = ((double)Aexc1 - (double)Aexc0)/(1.0/fexcDes[NfexcDes-1] - 1.0/fexcDes[0]);
altb 0:e2a7d7f91e49 120 this->bAexcDes = (double)Aexc0 - aAexcDes/fexcDes[0];
altb 0:e2a7d7f91e49 121
altb 0:e2a7d7f91e49 122 fnyq = 1.0/2.0/(double)Ts;
altb 0:e2a7d7f91e49 123 pi2 = 2.0*pi;
altb 0:e2a7d7f91e49 124 pi2Ts = pi2*(double)Ts;
altb 0:e2a7d7f91e49 125 piDiv2 = pi/2.0;
altb 0:e2a7d7f91e49 126
altb 0:e2a7d7f91e49 127 sU = (double*)malloc(3*sizeof(double));
altb 0:e2a7d7f91e49 128 sY = (double*)malloc(3*sizeof(double));
altb 0:e2a7d7f91e49 129 reset();
altb 0:e2a7d7f91e49 130 }
altb 0:e2a7d7f91e49 131
altb 0:e2a7d7f91e49 132 GPA::GPA(float f0, float f1, float *fexcDes, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1)
altb 0:e2a7d7f91e49 133 {
altb 0:e2a7d7f91e49 134 // convert fexcDes from float to double, it is assumed that it is sorted
altb 0:e2a7d7f91e49 135 this->NfexcDes = NfexcDes;
altb 0:e2a7d7f91e49 136 this->fexcDes = (double*)malloc(NfexcDes*sizeof(double));
altb 0:e2a7d7f91e49 137 for(int i = 0; i < NfexcDes; i++) {
altb 0:e2a7d7f91e49 138 this->fexcDes[i] = (double)fexcDes[i];
altb 0:e2a7d7f91e49 139 }
altb 0:e2a7d7f91e49 140 this->NperMin = NperMin;
altb 0:e2a7d7f91e49 141 this->NmeasMin = NmeasMin;
altb 0:e2a7d7f91e49 142 this->Ts = (double)Ts;
altb 0:e2a7d7f91e49 143
altb 0:e2a7d7f91e49 144 // calculate coefficients for decreasing amplitude (1/fexc)
altb 0:e2a7d7f91e49 145 this->aAexcDes = ((double)Aexc1 - (double)Aexc0)/(1.0/(double)f1 - 1.0/(double)f0);
altb 0:e2a7d7f91e49 146 this->bAexcDes = (double)Aexc0 - aAexcDes/fexcDes[0];
altb 0:e2a7d7f91e49 147
altb 0:e2a7d7f91e49 148 fnyq = 1.0/2.0/(double)Ts;
altb 0:e2a7d7f91e49 149 pi2 = 2.0*pi;
altb 0:e2a7d7f91e49 150 pi2Ts = pi2*(double)Ts;
altb 0:e2a7d7f91e49 151 piDiv2 = pi/2.0;
altb 0:e2a7d7f91e49 152
altb 0:e2a7d7f91e49 153 sU = (double*)malloc(3*sizeof(double));
altb 0:e2a7d7f91e49 154 sY = (double*)malloc(3*sizeof(double));
altb 0:e2a7d7f91e49 155 reset();
altb 0:e2a7d7f91e49 156 }
altb 0:e2a7d7f91e49 157
altb 0:e2a7d7f91e49 158 GPA::GPA(float *fexcDes, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1)
altb 0:e2a7d7f91e49 159 {
altb 0:e2a7d7f91e49 160 // convert fexcDes from float to double, it is assumed that it is sorted
altb 0:e2a7d7f91e49 161 this->NfexcDes = NfexcDes;
altb 0:e2a7d7f91e49 162 this->fexcDes = (double*)malloc(NfexcDes*sizeof(double));
altb 0:e2a7d7f91e49 163 for(int i = 0; i < NfexcDes; i++) {
altb 0:e2a7d7f91e49 164 this->fexcDes[i] = (double)fexcDes[i];
altb 0:e2a7d7f91e49 165 }
altb 0:e2a7d7f91e49 166 this->NperMin = NperMin;
altb 0:e2a7d7f91e49 167 this->NmeasMin = NmeasMin;
altb 0:e2a7d7f91e49 168 this->Ts = (double)Ts;
altb 0:e2a7d7f91e49 169
altb 0:e2a7d7f91e49 170 // calculate coefficients for decreasing amplitude (1/fexc)
altb 0:e2a7d7f91e49 171 this->aAexcDes = ((double)Aexc1 - (double)Aexc0)/(1.0/this->fexcDes[NfexcDes-1] - 1.0/this->fexcDes[0]);
altb 0:e2a7d7f91e49 172 this->bAexcDes = (double)Aexc0 - aAexcDes/fexcDes[0];
altb 0:e2a7d7f91e49 173
altb 0:e2a7d7f91e49 174 fnyq = 1.0/2.0/(double)Ts;
altb 0:e2a7d7f91e49 175 pi2 = 2.0*pi;
altb 0:e2a7d7f91e49 176 pi2Ts = pi2*(double)Ts;
altb 0:e2a7d7f91e49 177 piDiv2 = pi/2.0;
altb 0:e2a7d7f91e49 178
altb 0:e2a7d7f91e49 179 sU = (double*)malloc(3*sizeof(double));
altb 0:e2a7d7f91e49 180 sY = (double*)malloc(3*sizeof(double));
altb 0:e2a7d7f91e49 181 reset();
altb 0:e2a7d7f91e49 182 }
altb 0:e2a7d7f91e49 183
altb 0:e2a7d7f91e49 184 GPA::~GPA() {}
altb 0:e2a7d7f91e49 185
altb 0:e2a7d7f91e49 186 void GPA::reset()
altb 0:e2a7d7f91e49 187 {
altb 0:e2a7d7f91e49 188 Nmeas = 0;
altb 0:e2a7d7f91e49 189 Nper = 0;
altb 0:e2a7d7f91e49 190 fexc = 0.0;
altb 0:e2a7d7f91e49 191 fexcPast = 0.0;
altb 0:e2a7d7f91e49 192 ii = 1; // iterating through desired frequency points
altb 0:e2a7d7f91e49 193 jj = 1; // iterating through measurement points w.r.t. reachable frequency
altb 0:e2a7d7f91e49 194 scaleG = 0.0;
altb 0:e2a7d7f91e49 195 cr = 0.0;
altb 0:e2a7d7f91e49 196 ci = 0.0;
altb 0:e2a7d7f91e49 197 for(int i = 0; i < 3; i++) {
altb 0:e2a7d7f91e49 198 sU[i] = 0.0;
altb 0:e2a7d7f91e49 199 sY[i] = 0.0;
altb 0:e2a7d7f91e49 200 }
altb 0:e2a7d7f91e49 201 sinarg = 0.0;
altb 0:e2a7d7f91e49 202 NmeasTotal = 0;
altb 0:e2a7d7f91e49 203 Aexc = 0.0;
altb 0:e2a7d7f91e49 204 pi2Tsfexc = 0.0;
altb 0:e2a7d7f91e49 205 }
altb 0:e2a7d7f91e49 206
altb 0:e2a7d7f91e49 207 float GPA::update(double inp, double out)
altb 0:e2a7d7f91e49 208 {
altb 0:e2a7d7f91e49 209 // a new frequency point has been reached
altb 0:e2a7d7f91e49 210 if(jj == 1) {
altb 0:e2a7d7f91e49 211 // get a new unique frequency point
altb 0:e2a7d7f91e49 212 while(fexc == fexcPast) {
altb 0:e2a7d7f91e49 213 // measurement finished
altb 0:e2a7d7f91e49 214 if(ii > NfexcDes) {
altb 0:e2a7d7f91e49 215 return 0.0f;
altb 0:e2a7d7f91e49 216 }
altb 0:e2a7d7f91e49 217 calcGPAmeasPara(fexcDes[ii - 1]);
altb 0:e2a7d7f91e49 218 // secure fexc is not higher or equal to nyquist frequency
altb 0:e2a7d7f91e49 219 if(fexc >= fnyq) {
altb 0:e2a7d7f91e49 220 fexc = fexcPast;
altb 0:e2a7d7f91e49 221 }
altb 0:e2a7d7f91e49 222 // no frequency found
altb 0:e2a7d7f91e49 223 if(fexc == fexcPast) {
altb 0:e2a7d7f91e49 224 ii += 1;
altb 0:e2a7d7f91e49 225 } else {
altb 0:e2a7d7f91e49 226 Aexc = aAexcDes/fexc + bAexcDes;
altb 0:e2a7d7f91e49 227 pi2Tsfexc = pi2Ts*fexc;
altb 0:e2a7d7f91e49 228 }
altb 0:e2a7d7f91e49 229 }
altb 0:e2a7d7f91e49 230 // secure sinarg starts at 0 (numerically maybe not given)
altb 0:e2a7d7f91e49 231 sinarg = 0.0;
altb 0:e2a7d7f91e49 232 // filter scaling
altb 0:e2a7d7f91e49 233 scaleG = 1.0/sqrt((double)Nmeas);
altb 0:e2a7d7f91e49 234 // filter coefficients
altb 0:e2a7d7f91e49 235 cr = cos(pi2Tsfexc);
altb 0:e2a7d7f91e49 236 ci = sin(pi2Tsfexc);
altb 0:e2a7d7f91e49 237 // filter storage
altb 0:e2a7d7f91e49 238 for(int i = 0; i < 3; i++) {
altb 0:e2a7d7f91e49 239 sU[i] = 0.0;
altb 0:e2a7d7f91e49 240 sY[i] = 0.0;
altb 0:e2a7d7f91e49 241 }
altb 0:e2a7d7f91e49 242 }
altb 0:e2a7d7f91e49 243 // filter step for signal su
altb 0:e2a7d7f91e49 244 sU[0] = scaleG*inp + 2.0*cr*sU[1] - sU[2];
altb 0:e2a7d7f91e49 245 sU[2] = sU[1];
altb 0:e2a7d7f91e49 246 sU[1] = sU[0];
altb 0:e2a7d7f91e49 247 // filter step for signal sy
altb 0:e2a7d7f91e49 248 sY[0] = scaleG*out + 2.0*cr*sY[1] - sY[2];
altb 0:e2a7d7f91e49 249 sY[2] = sY[1];
altb 0:e2a7d7f91e49 250 sY[1] = sY[0];
altb 0:e2a7d7f91e49 251 // measurement of frequencypoint is finished
altb 0:e2a7d7f91e49 252 if(jj == Nmeas) {
altb 0:e2a7d7f91e49 253 jj = 1;
altb 0:e2a7d7f91e49 254 ii += 1;
altb 0:e2a7d7f91e49 255 fexcPast = fexc;
altb 0:e2a7d7f91e49 256 // calculate the one point dft
altb 0:e2a7d7f91e49 257 double Ureal = 2.0*scaleG*(cr*sU[1] - sU[2]);
altb 0:e2a7d7f91e49 258 double Uimag = 2.0*scaleG*ci*sU[1];
altb 0:e2a7d7f91e49 259 double Yreal = 2.0*scaleG*(cr*sY[1] - sY[2]);
altb 0:e2a7d7f91e49 260 double Yimag = 2.0*scaleG*ci*sY[1];
altb 0:e2a7d7f91e49 261 // calculate magnitude and angle
altb 0:e2a7d7f91e49 262 float Umag = (float)(sqrt(Ureal*Ureal + Uimag*Uimag));
altb 0:e2a7d7f91e49 263 float Ymag = (float)(sqrt(Yreal*Yreal + Yimag*Yimag));
altb 0:e2a7d7f91e49 264 float absGyu = (float)(Ymag/Umag);
altb 0:e2a7d7f91e49 265 float angGyu = (float)(atan2(Yimag, Yreal) - atan2(Uimag, Ureal));
altb 0:e2a7d7f91e49 266 float absGye = (float)(Ymag/Aexc);
altb 0:e2a7d7f91e49 267 float angGye = (float)(atan2(Yimag, Yreal) + piDiv2);
altb 0:e2a7d7f91e49 268 // user info
altb 0:e2a7d7f91e49 269 if(ii == 2) {
altb 0:e2a7d7f91e49 270 printLine();
altb 0:e2a7d7f91e49 271 printf(" fexc[Hz] |Gyu| ang(Gyu) |Gye| ang(Gye) |E| |U| |Y|\r\n");
altb 0:e2a7d7f91e49 272 printLine();
altb 0:e2a7d7f91e49 273 }
altb 0:e2a7d7f91e49 274 printf("%11.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\r\n", (float)fexc, absGyu, angGyu, absGye, angGye, (float)Aexc, Umag, Ymag);
altb 0:e2a7d7f91e49 275 } else {
altb 0:e2a7d7f91e49 276 jj += 1;
altb 0:e2a7d7f91e49 277 }
altb 0:e2a7d7f91e49 278 sinarg = fmod(sinarg + pi2Tsfexc, pi2);
altb 0:e2a7d7f91e49 279 NmeasTotal += 1;
altb 0:e2a7d7f91e49 280 return (float)(Aexc*sin(sinarg));
altb 0:e2a7d7f91e49 281 }
altb 0:e2a7d7f91e49 282
altb 0:e2a7d7f91e49 283 void GPA::fexcDesLogspace(double fMin, double fMax, int NfexcDes)
altb 0:e2a7d7f91e49 284 {
altb 0:e2a7d7f91e49 285 // calculate logarithmic spaced frequency points
altb 0:e2a7d7f91e49 286 double Gain = log10(fMax/fMin)/((double)NfexcDes - 1.0);
altb 0:e2a7d7f91e49 287 double expon = 0.0;
altb 0:e2a7d7f91e49 288 for(int i = 0; i < NfexcDes; i++) {
altb 0:e2a7d7f91e49 289 fexcDes[i] = fMin*pow(10.0, expon);
altb 0:e2a7d7f91e49 290 expon += Gain;
altb 0:e2a7d7f91e49 291 }
altb 0:e2a7d7f91e49 292 }
altb 0:e2a7d7f91e49 293
altb 0:e2a7d7f91e49 294 void GPA::calcGPAmeasPara(double fexcDes_i)
altb 0:e2a7d7f91e49 295 {
altb 0:e2a7d7f91e49 296 // Nmeas has to be an integer
altb 0:e2a7d7f91e49 297 Nper = NperMin;
altb 0:e2a7d7f91e49 298 Nmeas = (int)floor((double)Nper/fexcDes_i/Ts + 0.5);
altb 0:e2a7d7f91e49 299 // secure that the minimal number of measurements is fullfilled
altb 0:e2a7d7f91e49 300 int Ndelta = NmeasMin - Nmeas;
altb 0:e2a7d7f91e49 301 if(Ndelta > 0) {
altb 0:e2a7d7f91e49 302 Nper = (int)ceil((double)NmeasMin*fexcDes_i*Ts);
altb 0:e2a7d7f91e49 303 Nmeas = (int)floor((double)Nper/fexcDes_i/Ts + 0.5);
altb 0:e2a7d7f91e49 304 }
altb 0:e2a7d7f91e49 305 // evaluating reachable frequency
altb 0:e2a7d7f91e49 306 fexc = (double)Nper/(double)Nmeas/Ts;
altb 0:e2a7d7f91e49 307 }
altb 0:e2a7d7f91e49 308
altb 0:e2a7d7f91e49 309 void GPA::printLine()
altb 0:e2a7d7f91e49 310 {
altb 0:e2a7d7f91e49 311 printf("-----------------------------------------------------------------------------------------\r\n");
altb 0:e2a7d7f91e49 312 }
altb 0:e2a7d7f91e49 313
altb 0:e2a7d7f91e49 314 void GPA::printGPAfexcDes()
altb 0:e2a7d7f91e49 315 {
altb 0:e2a7d7f91e49 316 printLine();
altb 0:e2a7d7f91e49 317 for(int i = 0; i < NfexcDes; i++) {
altb 0:e2a7d7f91e49 318 printf("%9.4f\r\n", (float)fexcDes[i]);
altb 0:e2a7d7f91e49 319 }
altb 0:e2a7d7f91e49 320 }
altb 0:e2a7d7f91e49 321
altb 0:e2a7d7f91e49 322 void GPA::printGPAmeasPara()
altb 0:e2a7d7f91e49 323 {
altb 0:e2a7d7f91e49 324 printLine();
altb 0:e2a7d7f91e49 325 printf(" fexcDes[Hz] fexc[Hz] Aexc Nmeas Nper\r\n");
altb 0:e2a7d7f91e49 326 printLine();
altb 0:e2a7d7f91e49 327 for(int i = 0; i < NfexcDes; i++) {
altb 0:e2a7d7f91e49 328 calcGPAmeasPara(fexcDes[i]);
altb 0:e2a7d7f91e49 329 if(fexc == fexcPast || fexc >= fnyq) {
altb 0:e2a7d7f91e49 330 fexc = 0.0;
altb 0:e2a7d7f91e49 331 Nmeas = 0;
altb 0:e2a7d7f91e49 332 Nper = 0;
altb 0:e2a7d7f91e49 333 Aexc = 0.0;
altb 0:e2a7d7f91e49 334 } else {
altb 0:e2a7d7f91e49 335 Aexc = aAexcDes/fexc + bAexcDes;
altb 0:e2a7d7f91e49 336 fexcPast = fexc;
altb 0:e2a7d7f91e49 337 }
altb 0:e2a7d7f91e49 338 NmeasTotal += Nmeas;
altb 0:e2a7d7f91e49 339 printf("%12.2e %9.2e %10.2e %7i %6i \r\n", (float)fexcDes[i], (float)fexc, (float)Aexc, Nmeas, Nper);
altb 0:e2a7d7f91e49 340 }
altb 0:e2a7d7f91e49 341 printGPAmeasTime();
altb 0:e2a7d7f91e49 342 reset();
altb 0:e2a7d7f91e49 343 }
altb 0:e2a7d7f91e49 344
altb 0:e2a7d7f91e49 345 void GPA::printGPAmeasTime()
altb 0:e2a7d7f91e49 346 {
altb 0:e2a7d7f91e49 347 printLine();
altb 0:e2a7d7f91e49 348 printf(" number of data points: %9i\r\n", NmeasTotal);
altb 0:e2a7d7f91e49 349 printf(" measurment time in sec: %9.2f\r\n", (float)((double)NmeasTotal*Ts));
altb 0:e2a7d7f91e49 350 }
altb 0:e2a7d7f91e49 351
altb 0:e2a7d7f91e49 352 void GPA::printNfexcDes()
altb 0:e2a7d7f91e49 353 {
altb 0:e2a7d7f91e49 354 printLine();
altb 0:e2a7d7f91e49 355 printf(" number of frequancy points: %2i\r\n", NfexcDes);
altb 0:e2a7d7f91e49 356 }