Fertig
Dependencies: mbed
Fork of RT2_P3_students by
GPA.cpp@0:78ca29b4c49e, 2018-04-03 (annotated)
- Committer:
- altb
- Date:
- Tue Apr 03 08:47:41 2018 +0000
- Revision:
- 0:78ca29b4c49e
- Child:
- 3:769ce5f06d3e
Cuboid Lab RT2 FS 2018
Who changed what in which revision?
User | Revision | Line number | New contents of line |
---|---|---|---|
altb | 0:78ca29b4c49e | 1 | /* |
altb | 0:78ca29b4c49e | 2 | GPA: frequency point wise gain and phase analyser to measure the frequency |
altb | 0:78ca29b4c49e | 3 | respone of dynamical system |
altb | 0:78ca29b4c49e | 4 | |
altb | 0:78ca29b4c49e | 5 | hint: the measurements should only be perfomed in closed loop |
altb | 0:78ca29b4c49e | 6 | assumption: the system is at the desired steady state of interest when |
altb | 0:78ca29b4c49e | 7 | the measurment starts |
altb | 0:78ca29b4c49e | 8 | |
altb | 0:78ca29b4c49e | 9 | exc(2) C: controller |
altb | 0:78ca29b4c49e | 10 | | P: plant |
altb | 0:78ca29b4c49e | 11 | v |
altb | 0:78ca29b4c49e | 12 | exc(1) --> o ->| C |--->o------->| P |----------> out |
altb | 0:78ca29b4c49e | 13 | ^ | | |
altb | 0:78ca29b4c49e | 14 | | --> inp | exc: excitation signal (E) |
altb | 0:78ca29b4c49e | 15 | | | inp: input plant (U) |
altb | 0:78ca29b4c49e | 16 | -------------------------------- out: output plant (Y) |
altb | 0:78ca29b4c49e | 17 | |
altb | 0:78ca29b4c49e | 18 | instantiate option 1: |
altb | 0:78ca29b4c49e | 19 | GPA(float fMin, float fMax, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1) |
altb | 0:78ca29b4c49e | 20 | |
altb | 0:78ca29b4c49e | 21 | fMin: minimal desired frequency that should be measured in Hz |
altb | 0:78ca29b4c49e | 22 | fMax: maximal desired frequency that should be measured in Hz |
altb | 0:78ca29b4c49e | 23 | NfexcDes: number of logarithmic equaly spaced frequency points |
altb | 0:78ca29b4c49e | 24 | NperMin: minimal number of periods that are used for each frequency point |
altb | 0:78ca29b4c49e | 25 | NmeasMin: maximal number of samples that are used for each frequency point |
altb | 0:78ca29b4c49e | 26 | Ts: sampling time |
altb | 0:78ca29b4c49e | 27 | Aexc0: excitation amplitude at fMin |
altb | 0:78ca29b4c49e | 28 | Aexc1: excitation amplitude at fMax |
altb | 0:78ca29b4c49e | 29 | |
altb | 0:78ca29b4c49e | 30 | hints: the amplitude drops with 1/fexc, if you're using |
altb | 0:78ca29b4c49e | 31 | Axc1 = Aexc0/fMax the d/dt exc = const., this is recommended |
altb | 0:78ca29b4c49e | 32 | if your controller does not have a rolloff. |
altb | 0:78ca29b4c49e | 33 | if a desired frequency point is not measured try to increase Nmeas. |
altb | 0:78ca29b4c49e | 34 | |
altb | 0:78ca29b4c49e | 35 | pseudo code for a closed loop measurement with a proportional controller Kp: |
altb | 0:78ca29b4c49e | 36 | |
altb | 0:78ca29b4c49e | 37 | float inp = "measurement of inputsignal"; |
altb | 0:78ca29b4c49e | 38 | float out = "mesurement of outputsignal"; |
altb | 0:78ca29b4c49e | 39 | float exc = myGPA(inp, out); |
altb | 0:78ca29b4c49e | 40 | float off = "desired steady state off the system"; |
altb | 0:78ca29b4c49e | 41 | |
altb | 0:78ca29b4c49e | 42 | excitation input at (1): |
altb | 0:78ca29b4c49e | 43 | inp = Kp*(exc + off - out); |
altb | 0:78ca29b4c49e | 44 | |
altb | 0:78ca29b4c49e | 45 | excitation input at (2): |
altb | 0:78ca29b4c49e | 46 | inp = Kp*(off - out) + exc; |
altb | 0:78ca29b4c49e | 47 | |
altb | 0:78ca29b4c49e | 48 | usage: |
altb | 0:78ca29b4c49e | 49 | exc = myGPA(inp, out) does update the internal states of the gpa at the |
altb | 0:78ca29b4c49e | 50 | timestep k and returns the excitation signal for the timestep k+1. the |
altb | 0:78ca29b4c49e | 51 | results are plotted to a terminal (putty) over serial cennection and look |
altb | 0:78ca29b4c49e | 52 | as follows: |
altb | 0:78ca29b4c49e | 53 | ----------------------------------------------------------------------------------------- |
altb | 0:78ca29b4c49e | 54 | fexc[Hz] |Gyu| ang(Gyu) |Gye| ang(Gye) |E| |U| |Y| |
altb | 0:78ca29b4c49e | 55 | ----------------------------------------------------------------------------------------- |
altb | 0:78ca29b4c49e | 56 | 7.01e-01 1.08e+01 -7.45e-02 1.08e+01 -7.38e-02 9.99e-01 9.99e-01 1.08e+01 |
altb | 0:78ca29b4c49e | 57 | |
altb | 0:78ca29b4c49e | 58 | in matlab you can use: |
altb | 0:78ca29b4c49e | 59 | dataNucleo = [... insert measurement data here ...]; |
altb | 0:78ca29b4c49e | 60 | g = frd(dataNucleo(:,2).*exp(1i*dataNucleo(:,3)), dataNucleo(:,1), Ts, 'Units', 'Hz'); |
altb | 0:78ca29b4c49e | 61 | gcl = frd(dataNucleo(:,4).*exp(1i*dataNucleo(:,5)), dataNucleo(:,1), Ts, 'Units', 'Hz'); |
altb | 0:78ca29b4c49e | 62 | |
altb | 0:78ca29b4c49e | 63 | if you're evaluating more than one measurement which contain equal frequency points try: |
altb | 0:78ca29b4c49e | 64 | dataNucleo = [dataNucleo1; dataNucleo2]; |
altb | 0:78ca29b4c49e | 65 | [~, ind] = unique(dataNucleo(:,1),'stable'); |
altb | 0:78ca29b4c49e | 66 | dataNucleo = dataNucleo(ind,:); |
altb | 0:78ca29b4c49e | 67 | |
altb | 0:78ca29b4c49e | 68 | autor: M.E. Peter |
altb | 0:78ca29b4c49e | 69 | */ |
altb | 0:78ca29b4c49e | 70 | |
altb | 0:78ca29b4c49e | 71 | #include "GPA.h" |
altb | 0:78ca29b4c49e | 72 | #include "mbed.h" |
altb | 0:78ca29b4c49e | 73 | #include "math.h" |
altb | 0:78ca29b4c49e | 74 | #define pi 3.1415927f |
altb | 0:78ca29b4c49e | 75 | |
altb | 0:78ca29b4c49e | 76 | using namespace std; |
altb | 0:78ca29b4c49e | 77 | |
altb | 0:78ca29b4c49e | 78 | GPA::GPA(float fMin, float fMax, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1) |
altb | 0:78ca29b4c49e | 79 | { |
altb | 0:78ca29b4c49e | 80 | this->NfexcDes = NfexcDes; |
altb | 0:78ca29b4c49e | 81 | this->NperMin = NperMin; |
altb | 0:78ca29b4c49e | 82 | this->NmeasMin = NmeasMin; |
altb | 0:78ca29b4c49e | 83 | this->Ts = Ts; |
altb | 0:78ca29b4c49e | 84 | |
altb | 0:78ca29b4c49e | 85 | // calculate logarithmic spaced frequency points |
altb | 0:78ca29b4c49e | 86 | fexcDes = (float*)malloc(NfexcDes*sizeof(float)); |
altb | 0:78ca29b4c49e | 87 | fexcDesLogspace(fMin, fMax, NfexcDes); |
altb | 0:78ca29b4c49e | 88 | |
altb | 0:78ca29b4c49e | 89 | // calculate coefficients for decreasing amplitude (1/fexc) |
altb | 0:78ca29b4c49e | 90 | this->aAexcDes = (Aexc1 - Aexc0)/(1.0f/fexcDes[NfexcDes-1] - 1.0f/fexcDes[0]); |
altb | 0:78ca29b4c49e | 91 | this->bAexcDes = Aexc0 - aAexcDes/fexcDes[0]; |
altb | 0:78ca29b4c49e | 92 | |
altb | 0:78ca29b4c49e | 93 | fnyq = 1/2.0f/Ts; |
altb | 0:78ca29b4c49e | 94 | pi2 = 2.0f*pi; |
altb | 0:78ca29b4c49e | 95 | pi2Ts = pi2*Ts; |
altb | 0:78ca29b4c49e | 96 | piDiv2 = pi/2.0f; |
altb | 0:78ca29b4c49e | 97 | |
altb | 0:78ca29b4c49e | 98 | sU = (float*)malloc(3*sizeof(float)); |
altb | 0:78ca29b4c49e | 99 | sY = (float*)malloc(3*sizeof(float)); |
altb | 0:78ca29b4c49e | 100 | reset(); |
altb | 0:78ca29b4c49e | 101 | } |
altb | 0:78ca29b4c49e | 102 | |
altb | 0:78ca29b4c49e | 103 | GPA::~GPA() {} |
altb | 0:78ca29b4c49e | 104 | |
altb | 0:78ca29b4c49e | 105 | void GPA::reset() |
altb | 0:78ca29b4c49e | 106 | { |
altb | 0:78ca29b4c49e | 107 | Nmeas = 0; |
altb | 0:78ca29b4c49e | 108 | Nper = 0; |
altb | 0:78ca29b4c49e | 109 | fexc = 0.0f; |
altb | 0:78ca29b4c49e | 110 | fexcPast = 0.0f; |
altb | 0:78ca29b4c49e | 111 | ii = 1; // iterating through desired frequency points |
altb | 0:78ca29b4c49e | 112 | jj = 1; // iterating through measurement points w.r.t. reachable frequency |
altb | 0:78ca29b4c49e | 113 | scaleG = 0.0f; |
altb | 0:78ca29b4c49e | 114 | cr = 0.0f; |
altb | 0:78ca29b4c49e | 115 | ci = 0.0f; |
altb | 0:78ca29b4c49e | 116 | for(int i = 0; i < 3; i++) { |
altb | 0:78ca29b4c49e | 117 | sU[i] = 0.0f; |
altb | 0:78ca29b4c49e | 118 | sY[i] = 0.0f; |
altb | 0:78ca29b4c49e | 119 | } |
altb | 0:78ca29b4c49e | 120 | sinarg = 0.0f; |
altb | 0:78ca29b4c49e | 121 | NmeasTotal = 0; |
altb | 0:78ca29b4c49e | 122 | Aexc = 0.0f; |
altb | 0:78ca29b4c49e | 123 | pi2Tsfexc = 0.0f; |
altb | 0:78ca29b4c49e | 124 | } |
altb | 0:78ca29b4c49e | 125 | |
altb | 0:78ca29b4c49e | 126 | float GPA::update(float inp, float out) |
altb | 0:78ca29b4c49e | 127 | { |
altb | 0:78ca29b4c49e | 128 | // a new frequency point has been reached |
altb | 0:78ca29b4c49e | 129 | if(jj == 1) { |
altb | 0:78ca29b4c49e | 130 | // get a new unique frequency point |
altb | 0:78ca29b4c49e | 131 | while(fexc == fexcPast) { |
altb | 0:78ca29b4c49e | 132 | // measurement finished |
altb | 0:78ca29b4c49e | 133 | if(ii > NfexcDes) { |
altb | 0:78ca29b4c49e | 134 | return 0.0f; |
altb | 0:78ca29b4c49e | 135 | } |
altb | 0:78ca29b4c49e | 136 | calcGPAmeasPara(fexcDes[ii - 1]); |
altb | 0:78ca29b4c49e | 137 | // secure fexc is not higher or equal to nyquist frequency |
altb | 0:78ca29b4c49e | 138 | if(fexc >= fnyq) { |
altb | 0:78ca29b4c49e | 139 | fexc = fexcPast; |
altb | 0:78ca29b4c49e | 140 | } |
altb | 0:78ca29b4c49e | 141 | // no frequency found |
altb | 0:78ca29b4c49e | 142 | if(fexc == fexcPast) { |
altb | 0:78ca29b4c49e | 143 | ii += 1; |
altb | 0:78ca29b4c49e | 144 | } else { |
altb | 0:78ca29b4c49e | 145 | Aexc = aAexcDes/fexc + bAexcDes; |
altb | 0:78ca29b4c49e | 146 | pi2Tsfexc = pi2Ts*fexc; |
altb | 0:78ca29b4c49e | 147 | } |
altb | 0:78ca29b4c49e | 148 | } |
altb | 0:78ca29b4c49e | 149 | fexcPast = fexc; |
altb | 0:78ca29b4c49e | 150 | // filter scaling |
altb | 0:78ca29b4c49e | 151 | scaleG = 1.0f/sqrt((float)Nmeas); |
altb | 0:78ca29b4c49e | 152 | // filter coefficients |
altb | 0:78ca29b4c49e | 153 | cr = cos(pi2Tsfexc); |
altb | 0:78ca29b4c49e | 154 | ci = sin(pi2Tsfexc); |
altb | 0:78ca29b4c49e | 155 | // filter storage |
altb | 0:78ca29b4c49e | 156 | for(int i = 0; i < 3; i++) { |
altb | 0:78ca29b4c49e | 157 | sU[i] = 0.0f; |
altb | 0:78ca29b4c49e | 158 | sY[i] = 0.0f; |
altb | 0:78ca29b4c49e | 159 | } |
altb | 0:78ca29b4c49e | 160 | } |
altb | 0:78ca29b4c49e | 161 | // filter step for signal su |
altb | 0:78ca29b4c49e | 162 | sU[0] = scaleG*inp + 2.0f*cr*sU[1] - sU[2]; |
altb | 0:78ca29b4c49e | 163 | sU[2] = sU[1]; |
altb | 0:78ca29b4c49e | 164 | sU[1] = sU[0]; |
altb | 0:78ca29b4c49e | 165 | // filter step for signal sy |
altb | 0:78ca29b4c49e | 166 | sY[0] = scaleG*out + 2.0f*cr*sY[1] - sY[2]; |
altb | 0:78ca29b4c49e | 167 | sY[2] = sY[1]; |
altb | 0:78ca29b4c49e | 168 | sY[1] = sY[0]; |
altb | 0:78ca29b4c49e | 169 | // measurement of frequencypoint is finished |
altb | 0:78ca29b4c49e | 170 | if(jj == Nmeas) { |
altb | 0:78ca29b4c49e | 171 | jj = 1; |
altb | 0:78ca29b4c49e | 172 | ii += 1; |
altb | 0:78ca29b4c49e | 173 | // calculate the one point dft |
altb | 0:78ca29b4c49e | 174 | float Ureal = 2.0f*scaleG*(cr*sU[1] - sU[2]); |
altb | 0:78ca29b4c49e | 175 | float Uimag = 2.0f*scaleG*ci*sU[1]; |
altb | 0:78ca29b4c49e | 176 | float Yreal = 2.0f*scaleG*(cr*sY[1] - sY[2]); |
altb | 0:78ca29b4c49e | 177 | float Yimag = 2.0f*scaleG*ci*sY[1]; |
altb | 0:78ca29b4c49e | 178 | // calculate magnitude and angle |
altb | 0:78ca29b4c49e | 179 | float Umag = sqrt(Ureal*Ureal + Uimag*Uimag); |
altb | 0:78ca29b4c49e | 180 | float Ymag = sqrt(Yreal*Yreal + Yimag*Yimag); |
altb | 0:78ca29b4c49e | 181 | float absGyu = Ymag/Umag; |
altb | 0:78ca29b4c49e | 182 | float angGyu = atan2(Yimag, Yreal) - atan2(Uimag, Ureal); |
altb | 0:78ca29b4c49e | 183 | float absGye = Ymag/Aexc; |
altb | 0:78ca29b4c49e | 184 | float angGye = (atan2(Yimag, Yreal) + piDiv2); |
altb | 0:78ca29b4c49e | 185 | // user info |
altb | 0:78ca29b4c49e | 186 | if(ii == 1) { |
altb | 0:78ca29b4c49e | 187 | printLine(); |
altb | 0:78ca29b4c49e | 188 | printf(" fexc[Hz] |Gyu| ang(Gyu) |Gye| ang(Gye) |E| |U| |Y|\r\n"); |
altb | 0:78ca29b4c49e | 189 | printLine(); |
altb | 0:78ca29b4c49e | 190 | } |
altb | 0:78ca29b4c49e | 191 | printf("%11.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e %10.2e\r\n", fexc, absGyu, angGyu, absGye, angGye, Aexc, Umag, Ymag); |
altb | 0:78ca29b4c49e | 192 | } else { |
altb | 0:78ca29b4c49e | 193 | jj += 1; |
altb | 0:78ca29b4c49e | 194 | } |
altb | 0:78ca29b4c49e | 195 | sinarg = fmod(sinarg + pi2Tsfexc, pi2); |
altb | 0:78ca29b4c49e | 196 | NmeasTotal += 1; |
altb | 0:78ca29b4c49e | 197 | return Aexc*sin(sinarg); |
altb | 0:78ca29b4c49e | 198 | } |
altb | 0:78ca29b4c49e | 199 | |
altb | 0:78ca29b4c49e | 200 | void GPA::fexcDesLogspace(float fMin, float fMax, int NfexcDes) |
altb | 0:78ca29b4c49e | 201 | { |
altb | 0:78ca29b4c49e | 202 | // calculate logarithmic spaced frequency points |
altb | 0:78ca29b4c49e | 203 | float Gain = log10(fMax/fMin)/((float)NfexcDes - 1.0f); |
altb | 0:78ca29b4c49e | 204 | float expon = 0.0f; |
altb | 0:78ca29b4c49e | 205 | for(int i = 0; i < NfexcDes; i++) { |
altb | 0:78ca29b4c49e | 206 | fexcDes[i] = fMin*pow(10.0f, expon); |
altb | 0:78ca29b4c49e | 207 | expon += Gain; |
altb | 0:78ca29b4c49e | 208 | } |
altb | 0:78ca29b4c49e | 209 | } |
altb | 0:78ca29b4c49e | 210 | |
altb | 0:78ca29b4c49e | 211 | void GPA::calcGPAmeasPara(float fexcDes_i) |
altb | 0:78ca29b4c49e | 212 | { |
altb | 0:78ca29b4c49e | 213 | // Nmeas has to be an integer |
altb | 0:78ca29b4c49e | 214 | Nper = NperMin; |
altb | 0:78ca29b4c49e | 215 | Nmeas = (int)floor((float)Nper/fexcDes_i/Ts + 0.5f); |
altb | 0:78ca29b4c49e | 216 | // secure that the minimal number of measurements is fullfilled |
altb | 0:78ca29b4c49e | 217 | int Ndelta = NmeasMin - Nmeas; |
altb | 0:78ca29b4c49e | 218 | if(Ndelta > 0) { |
altb | 0:78ca29b4c49e | 219 | Nper = (int)ceil((float)NmeasMin*fexcDes_i*Ts); |
altb | 0:78ca29b4c49e | 220 | Nmeas = (int)floor((float)Nper/fexcDes_i/Ts + 0.5f); |
altb | 0:78ca29b4c49e | 221 | } |
altb | 0:78ca29b4c49e | 222 | // evaluating reachable frequency |
altb | 0:78ca29b4c49e | 223 | fexc = (float)Nper/(float)Nmeas/Ts; |
altb | 0:78ca29b4c49e | 224 | } |
altb | 0:78ca29b4c49e | 225 | |
altb | 0:78ca29b4c49e | 226 | void GPA::printLine() |
altb | 0:78ca29b4c49e | 227 | { |
altb | 0:78ca29b4c49e | 228 | printf("-----------------------------------------------------------------------------------------\r\n"); |
altb | 0:78ca29b4c49e | 229 | } |
altb | 0:78ca29b4c49e | 230 | |
altb | 0:78ca29b4c49e | 231 | void GPA::printGPAfexcDes() |
altb | 0:78ca29b4c49e | 232 | { |
altb | 0:78ca29b4c49e | 233 | printLine(); |
altb | 0:78ca29b4c49e | 234 | for(int i = 0; i < NfexcDes; i++) { |
altb | 0:78ca29b4c49e | 235 | printf("%9.4f\r\n", fexcDes[i]); |
altb | 0:78ca29b4c49e | 236 | } |
altb | 0:78ca29b4c49e | 237 | } |
altb | 0:78ca29b4c49e | 238 | |
altb | 0:78ca29b4c49e | 239 | void GPA::printGPAmeasPara() |
altb | 0:78ca29b4c49e | 240 | { |
altb | 0:78ca29b4c49e | 241 | printLine(); |
altb | 0:78ca29b4c49e | 242 | printf(" fexcDes[Hz] fexc[Hz] Aexc Nmeas Nper\r\n"); |
altb | 0:78ca29b4c49e | 243 | printLine(); |
altb | 0:78ca29b4c49e | 244 | for(int i = 0; i < NfexcDes; i++) { |
altb | 0:78ca29b4c49e | 245 | calcGPAmeasPara(fexcDes[i]); |
altb | 0:78ca29b4c49e | 246 | if(fexc == fexcPast || fexc >= fnyq) { |
altb | 0:78ca29b4c49e | 247 | fexc = 0.0f; |
altb | 0:78ca29b4c49e | 248 | Nmeas = 0; |
altb | 0:78ca29b4c49e | 249 | Nper = 0; |
altb | 0:78ca29b4c49e | 250 | Aexc = 0; |
altb | 0:78ca29b4c49e | 251 | } else { |
altb | 0:78ca29b4c49e | 252 | Aexc = aAexcDes/fexc + bAexcDes; |
altb | 0:78ca29b4c49e | 253 | fexcPast = fexc; |
altb | 0:78ca29b4c49e | 254 | } |
altb | 0:78ca29b4c49e | 255 | NmeasTotal += Nmeas; |
altb | 0:78ca29b4c49e | 256 | printf("%12.2e %9.2e %10.2e %7i %6i \r\n", fexcDes[i], fexc, Aexc, Nmeas, Nper); |
altb | 0:78ca29b4c49e | 257 | } |
altb | 0:78ca29b4c49e | 258 | printGPAmeasTime(); |
altb | 0:78ca29b4c49e | 259 | reset(); |
altb | 0:78ca29b4c49e | 260 | } |
altb | 0:78ca29b4c49e | 261 | |
altb | 0:78ca29b4c49e | 262 | void GPA::printGPAmeasTime() |
altb | 0:78ca29b4c49e | 263 | { |
altb | 0:78ca29b4c49e | 264 | printLine(); |
altb | 0:78ca29b4c49e | 265 | printf(" number of data points: %9i\r\n", NmeasTotal); |
altb | 0:78ca29b4c49e | 266 | printf(" measurment time in sec: %9.2f\r\n", (float)NmeasTotal*Ts); |
altb | 0:78ca29b4c49e | 267 | } |