Test program with the RT black boxes

Dependencies:   mbed

Embed: (wiki syntax)

« Back to documentation index

Show/hide line numbers GPA.cpp Source File

GPA.cpp

00001 /*
00002     GPA: Frequency point wise gain and phase analyser to measure the frequency respone function (FRF) of a dynamical system, based on the one point DFT
00003 
00004          Hint:        If the plant has a pole at zero, is unstable or weakly damped the measurement has to be perfomed 
00005                       in closed loop (this is NOT tfestimate, the algorithm is based on the one point DFT).
00006          Assumption:  The system is and remains at the desired steady state of interest when the measurment starts
00007 
00008     Instantiate option 0: ("Not a Jedi yet" users, for logarithmic equaly spaced frequency points)
00009     
00010         GPA(float fMin, float fMax, int NfexcDes, float Aexc0, float Aexc1, float Ts)
00011         
00012             fMin:       Minimal desired frequency that should be measured in Hz
00013             fMax:       Maximal desired frequency that should be measured in Hz
00014             NfexcDes:   Number of logarithmic equaly spaced frequency points between fMin and fMax
00015             Aexc0:      Excitation amplitude at fMin
00016             Aexc1:      Excitation amplitude at fMax
00017             Ts:         Sampling time in sec
00018             
00019             Default values are as follows:
00020             int NperMin   = 3;
00021             int NmeasMin  = (int)ceil(1.0f/Ts);
00022             int NstartMin = (int)ceil(3.0f/Ts);
00023             int NsweepMin = 0;
00024 
00025     Instantiate option 1: ("Jedi or Sith Lord", for logarithmic equaly spaced frequency points)
00026     
00027         GPA(float fMin, float fMax, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1, int NstartMin, int NsweepMin)
00028         
00029             fMin:       Minimal desired frequency that should be measured in Hz
00030             fMax:       Maximal desired frequency that should be measured in Hz
00031             NfexcDes:   Number of logarithmic equaly spaced frequency points
00032             NperMin:    Minimal number of periods that are used for each frequency point
00033             NmeasMin:   Minimal number of samples that are used for each frequency point
00034             Ts:         Sampling time in sec
00035             Aexc0:      Excitation amplitude at fMin
00036             Aexc1:      Excitation amplitude at fMax
00037             NstartMin:  Minimal number of samples to sweep to the first frequency point (can be equal 0)
00038             NsweepMin:  Minimal number of samples to sweep from freq. point to freq. point (can be equal 0)
00039 
00040  
00041     Instantiate option 2: (for a second, refined frequency grid measurement)
00042     
00043         GPA(float f0, float f1, float *fexcDes, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1, int NstartMin, int NsweepMin)
00044         
00045             f0:         Frequency point for the calculation of Aexc0 in Hz (should be chosen like in the first measurement)
00046             f1:         Frequency point for the calculation of Aexc1 in Hz (should be chosen like in the first measurement)
00047             *fexcDes:   Sorted frequency point array in Hz
00048             NfexcDes:   Length of fexcDes
00049             
00050             For the other parameters see above.
00051 
00052     Instantiate option 3: (for an arbitary but sorted frequency grid measurement)
00053     
00054         GPA(float *fexcDes, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1, int NstartMin, int NsweepMin)
00055         
00056             *fexcDes:   Sorted frequency point array in Hz
00057             Aexc0:      Excitation amplitude at fexcDes[0]
00058             Aexc1:      Excitation amplitude at fexcDes[NfexcDes-1]
00059             NfexcDes:   Length of fexcDes
00060             
00061             For the other parameters see above.
00062 
00063     Note:   The amplitude drops with 1/fexc, if you're using Axc1 = Aexc0/fMax then d/dt exc = const.,
00064             this is recommended if your controller does not have a rolloff. If a desired frequency point
00065             is not measured (could not be reached) try to increase Nmeas.
00066               
00067         
00068     Block diagram:
00069 
00070                 w (const.)    exc(2)                C: controller
00071                 |              |                    P: plant
00072                 v   e          v
00073      exc(1) --> o   ->| C |--->o------->| P |----------> out (y)
00074                 ^ -                |             |
00075                 |                   --> inp (u)  |  exc (R): excitation signal
00076                 |                                |  inp (U): input plant
00077                  --------------------------------   out (Y): output plant
00078               
00079               
00080     Pseudo code for an open loop measurement:
00081 
00082         - Measuring the plant P = Gyu = Gyr:
00083         
00084             u = w + exc;
00085             ... write output u here! it follows exc(k+1) ...
00086             exc = Wobble(exc, y);
00087             
00088             Closed loop FRF calculus with a stabilizing controller C:
00089                 S  = 1/(1 + C*P);  % ( exc1 -> e ,   1/(1 + C*P) ) sensitivity, contr. error rejection, robustness (1/max|S|)
00090                 T  = 1 - S;        % ( w -> y    , C*P/(1 + C*P) ) compl. sensitivity, tracking performance
00091                 CS = C*S;          % ( exc1 -> u ,   C/(1 + C*P) ) control effort, disturbance plant output on plant input
00092                 PS = P*S;          % ( exc2 -> y ,   P/(1 + C*P) ) compliance, disturbance plant input on plant output
00093 
00094 
00095     Pseudo code for a closed loop measurement with stabilizing controller C:
00096 
00097         Excitation at excitation input (1):   
00098         
00099         - Measuring the plant P = Gyu and the closed loop tf T = PC/(1 + PC) = Gyr:
00100         
00101             u = C(w - y + exc);
00102             ... write output u here! it follows exc(k+1) ...
00103             exc = Wobble(u, y);
00104                         
00105             Closed loop FRF calculus:
00106                 S  = 1 - T;
00107                 PS = P*S;
00108                 CS = T/P;
00109                 C  = C/S;
00110 
00111         Excitation at excitation input (2):
00112         
00113         - Measuring the plant P = Gyu and the closed loop tf PS = P/(1 + PC) = Gyr:
00114         
00115             u = C(w - y) + exc;
00116             ... write output u here! it follows exc(k+1) ...
00117             exc = Wobble(u, y);
00118         
00119             Closed loop FRF calculus:
00120                 S  = PS/P;
00121                 T  = 1 - S;
00122                 CS = T/P;
00123                 C  = C/S;
00124 
00125 
00126     Usage:
00127         exc(k+1) = myGPA(inp(k), out(k)) does update the internal states of the 
00128         gpa at the timestep k and returns the excitation signal for the timestep
00129         k+1. The FRF data are plotted to a terminal (Putty) over a serial 
00130         connection and look as follows:
00131         
00132 --------------------------------------------------------------------------------
00133   fexc[Hz]    |Gyu|    deg(Gyu)  |Gyr|    deg(Gyr)   |U|       |Y|       |R|
00134 --------------------------------------------------------------------------------
00135  5.0000e-02 1.001e+00   -0.309 1.001e+00   -0.309 4.000e-01 4.000e-01 4.005e-01
00136     .           .         .        .         .        .         .         .
00137     .           .         .        .         .        .         .         .
00138     .           .         .        .         .        .         .         .
00139   
00140     In Matlab you can use the editor as follows:
00141         data = [... insert measurement data here ...];
00142         gyu = frd(data(:,2).*exp(1i*data(:,3)*pi/180), data(:,1), Ts, 'Units', 'Hz');
00143         gyr = frd(data(:,4).*exp(1i*data(:,5)*pi/180), data(:,1), Ts, 'Units', 'Hz');
00144 
00145     If you're evaluating more than one measurement which contain equal frequency points use:
00146         data = [data1; data2];
00147         [~, ind] = unique(data(:,1), 'stable');
00148         data = data(ind,:);
00149 
00150 
00151     Autor and Copyrigth: 2018 / 2019 / M.E. Peter
00152     
00153 */
00154 
00155 #include "GPA.h"
00156 #include "mbed.h"
00157 #include "math.h"
00158 #define   pi 3.141592653589793
00159 
00160 using namespace std;
00161 
00162 // -----------------------------------------------------------------------------
00163 //      instantiate
00164 // -----------------------------------------------------------------------------
00165 
00166 GPA::GPA(float fMin, float fMax, int NfexcDes, float Aexc0, float Aexc1, float Ts)
00167 {
00168     int NperMin = 3;
00169     int NmeasMin = (int)ceil(1.0f/Ts);
00170     int NstartMin = (int)ceil(3.0f/Ts);
00171     int NsweepMin = 0;
00172     
00173     doPrint = true;
00174     assignParameters(NfexcDes, NperMin, NmeasMin, (double)Ts, NstartMin, NsweepMin);
00175 
00176     // calculate logarithmic spaced frequency points
00177     fexcDes = (double*)malloc(NfexcDes*sizeof(double));
00178     fexcDesLogspace((double)fMin, (double)fMax, NfexcDes);
00179 
00180     calculateDecreasingAmplitudeCoefficients((double)Aexc0, (double)Aexc1);
00181     initializeConstants((double)Ts);
00182     assignFilterStorage();
00183     reset();
00184 }
00185 
00186 GPA::GPA(float fMin, float fMax, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1, int NstartMin, int NsweepMin)
00187 {
00188     doPrint = true;
00189     assignParameters(NfexcDes, NperMin, NmeasMin, (double)Ts, NstartMin, NsweepMin);
00190 
00191     // calculate logarithmic spaced frequency points
00192     fexcDes = (double*)malloc(NfexcDes*sizeof(double));
00193     fexcDesLogspace((double)fMin, (double)fMax, NfexcDes);
00194 
00195     calculateDecreasingAmplitudeCoefficients((double)Aexc0, (double)Aexc1);
00196     initializeConstants((double)Ts);
00197     assignFilterStorage();
00198     reset();
00199 }
00200 
00201 GPA::GPA(float f0, float f1, float *fexcDes, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1, int NstartMin, int NsweepMin)
00202 {
00203     doPrint = true;
00204     assignParameters(NfexcDes, NperMin, NmeasMin, (double)Ts, NstartMin, NsweepMin);
00205     
00206     // convert fexcDes from float to double, it is assumed that it is sorted
00207     this->fexcDes = (double*)malloc(NfexcDes*sizeof(double));
00208     for(int i = 0; i < NfexcDes; i++) {
00209         this->fexcDes[i] = (double)fexcDes[i];
00210     }
00211     
00212     calculateDecreasingAmplitudeCoefficients((double)Aexc0, (double)Aexc1);
00213     initializeConstants((double)Ts);
00214     assignFilterStorage();
00215     reset();
00216 }
00217 
00218 GPA::GPA(float *fexcDes, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1, int NstartMin, int NsweepMin)
00219 {
00220     doPrint = true;
00221     assignParameters(NfexcDes, NperMin, NmeasMin, (double)Ts, NstartMin, NsweepMin);
00222     
00223     // convert fexcDes from float to double, it is assumed that it is sorted
00224     this->fexcDes = (double*)malloc(NfexcDes*sizeof(double));
00225     for(int i = 0; i < NfexcDes; i++) {
00226         this->fexcDes[i] = (double)fexcDes[i];
00227     }
00228 
00229     calculateDecreasingAmplitudeCoefficients((double)Aexc0, (double)Aexc1);
00230     initializeConstants((double)Ts);
00231     assignFilterStorage();
00232     reset();
00233 }
00234 
00235 GPA::GPA(float fMin, float fMax, int NfexcDes, int NperMin, int NmeasMin, float Ts, float Aexc0, float Aexc1, int NstartMin, int NsweepMin, bool doPrint)
00236 {
00237     this->doPrint = doPrint;
00238     assignParameters(NfexcDes, NperMin, NmeasMin, (double)Ts, NstartMin, NsweepMin);
00239 
00240     // calculate logarithmic spaced frequency points
00241     fexcDes = (double*)malloc(NfexcDes*sizeof(double));
00242     fexcDesLogspace((double)fMin, (double)fMax, NfexcDes);
00243 
00244     calculateDecreasingAmplitudeCoefficients((double)Aexc0, (double)Aexc1);
00245     initializeConstants((double)Ts);
00246     assignFilterStorage();
00247     reset();
00248 }
00249 
00250 // -----------------------------------------------------------------------------
00251 //      virtual and reset
00252 // -----------------------------------------------------------------------------
00253 
00254 GPA::~GPA() {}
00255 
00256 void GPA::reset()
00257 {
00258     Nmeas = 0;
00259     Nper = 0;
00260     dfexc = 0.0;
00261     fexc = 0.0;
00262     fexcPast = 0.0;
00263     i = 1; // iterating through desired frequency points
00264     j = 1; // iterating through measurement points w.r.t. reachable frequency
00265     scaleG = 0.0;
00266     cr = 0.0;
00267     ci = 0.0;
00268     for(int i = 0; i < 3; i++) {
00269         sU[i] = 0.0;
00270         sY[i] = 0.0;
00271     }
00272     sinarg = 0.0;
00273     NmeasTotal = 0;
00274     Aexc = 0.0;
00275     pi2Tsfexc = 0.0;
00276     Nsweep = NstartMin;
00277     bfexc = 0.0;
00278     afexc = 0.0;
00279     aAexc = 0.0;
00280     bAexc = 0.0;
00281     AexcOut = 0.0;
00282     gpaData.fexc = 0.0f;
00283     gpaData.absGyu = 0.0f;
00284     gpaData.angGyu = 0.0f;
00285     gpaData.absGyr = 0.0f;
00286     gpaData.angGyr = 0.0f;
00287     gpaData.Umag = 0.0f;
00288     gpaData.Ymag = 0.0f;
00289     gpaData.Rmag = 0.0f;
00290     gpaData.MeasPointFinished = false;
00291     
00292 }
00293 
00294 // -----------------------------------------------------------------------------
00295 //      update (operator)
00296 // -----------------------------------------------------------------------------  
00297 
00298 float GPA::update(double inp, double out)
00299 {
00300     // a new frequency point has been reached
00301     if(j == 1) {
00302         // user info
00303         if(i == 1 & doPrint) {
00304             printLine();
00305             printf("  fexc[Hz]    |Gyu|    deg(Gyu)  |Gyr|    deg(Gyr)   |U|       |Y|       |R|\r\n");
00306             printLine();
00307         }
00308         // get a new unique frequency point
00309         while(fexc == fexcPast) {
00310             // measurement finished
00311             if(i > NfexcDes) {
00312                 gpaData.MeasPointFinished = false;
00313                 return 0.0f;
00314             }
00315             calcGPAmeasPara(fexcDes[i - 1]);
00316             // secure fexc is not higher or equal to nyquist frequency
00317             if(fexc >= fnyq) {
00318                 fexc = fexcPast;
00319             }
00320             // no frequency found
00321             if(fexc == fexcPast) {
00322                 i += 1;
00323             } else {
00324                 Aexc = aAexcDes/fexc + bAexcDes;
00325                 pi2Tsfexc = pi2Ts*fexc;
00326             }
00327         }
00328         // filter scaling
00329         scaleG = 1.0/sqrt((double)Nmeas);
00330         // filter coefficients
00331         cr = cos(pi2Tsfexc);
00332         ci = sin(pi2Tsfexc);
00333         // set filter storage zero
00334         for(int i = 0; i < 3; i++) {
00335             sU[i] = 0.0;
00336             sY[i] = 0.0;
00337         }
00338         // calculate the parameters for the frequency sweep from fexcPast to fexc
00339         if(Nsweep > 0) calcGPAsweepPara();
00340         gpaData.MeasPointFinished = false;
00341     }
00342     // perfomre the sweep or measure
00343     if(j <= Nsweep) {
00344         dfexc = afexc*(double)j + bfexc;
00345         AexcOut = aAexc*(double)j + bAexc;   
00346     } else {
00347         dfexc = fexc;
00348         AexcOut = Aexc;
00349         // one point DFT filter step for signal su
00350         sU[0] = scaleG*inp + 2.0*cr*sU[1] - sU[2];
00351         sU[2] = sU[1];
00352         sU[1] = sU[0];
00353         // one point DFT filter step for signal sy
00354         sY[0] = scaleG*out + 2.0*cr*sY[1] - sY[2];
00355         sY[2] = sY[1];
00356         sY[1] = sY[0];
00357     }
00358     // secure sinarg starts at 0 (numerically maybe not given)
00359     if(j == 1 || j == Nsweep + 1) sinarg = 0.0;
00360     // measurement of frequencypoint is finished
00361     if(j == Nmeas + Nsweep) {
00362         fexcPast = fexc;
00363         AexcPast = Aexc;
00364         Nsweep = NsweepMin;
00365         // calculate the one point dft
00366         double Ureal = 2.0*scaleG*(cr*sU[1] - sU[2]);
00367         double Uimag = 2.0*scaleG*ci*sU[1];
00368         double Yreal = 2.0*scaleG*(cr*sY[1] - sY[2]);
00369         double Yimag = 2.0*scaleG*ci*sY[1];
00370         // calculate magnitude and angle
00371         double Umag = sqrt(Ureal*Ureal + Uimag*Uimag);
00372         double Ymag = sqrt(Yreal*Yreal + Yimag*Yimag);
00373         gpaData.fexc = (float)fexc;
00374         gpaData.absGyu = (float)(Ymag/Umag);
00375         gpaData.angGyu = (float)(wrapAngle(atan2(Yimag, Yreal) - atan2(Uimag, Ureal))*rad2deg);
00376         gpaData.absGyr = (float)(Ymag/Aexc);
00377         gpaData.angGyr = (float)(wrapAngle(atan2(Yimag, Yreal) + piDiv2)*rad2deg);
00378         gpaData.Umag = (float)(sqrt(Ureal*Ureal + Uimag*Uimag));
00379         gpaData.Ymag = (float)(sqrt(Yreal*Yreal + Yimag*Yimag));
00380         gpaData.Rmag = (float)Aexc;
00381         gpaData.MeasPointFinished = true;
00382         // user info
00383         if(doPrint) {
00384         printf("%11.4e %9.3e %8.3f %9.3e %8.3f %9.3e %9.3e %9.3e\r\n", gpaData.fexc, gpaData.absGyu, gpaData.angGyu, gpaData.absGyr, gpaData.angGyr, gpaData.Umag, gpaData.Ymag, gpaData.Rmag);
00385         }
00386         i += 1;
00387         j = 1;
00388     } else {
00389         j += 1;
00390     }
00391     // calculate the excitation
00392     sinarg = fmod(sinarg + pi2Ts*dfexc, pi2);
00393     NmeasTotal += 1;
00394     return (float)(AexcOut*sin(sinarg));
00395 }
00396 
00397 // -----------------------------------------------------------------------------
00398 //      private functions
00399 // -----------------------------------------------------------------------------  
00400 
00401 void GPA::assignParameters(int NfexcDes, int NperMin, int NmeasMin, double Ts, int NstartMin, int NsweepMin)
00402 {
00403     this->NfexcDes = NfexcDes;
00404     this->NperMin = NperMin;
00405     this->NmeasMin = NmeasMin;
00406     this->Ts = Ts;
00407     this->NstartMin = NstartMin;
00408     this->NsweepMin = NsweepMin;   
00409 }
00410 
00411 void GPA::calculateDecreasingAmplitudeCoefficients(double Aexc0, double Aexc1)
00412 {
00413     // calculate coefficients for decreasing amplitude (1/fexc)
00414     this->aAexcDes = (Aexc1 - Aexc0)/(1.0/fexcDes[NfexcDes-1] - 1.0/fexcDes[0]);
00415     this->bAexcDes = Aexc0 - aAexcDes/fexcDes[0];   
00416 }
00417 
00418 void GPA::initializeConstants(double Ts)
00419 {
00420     fnyq = 1.0/2.0/Ts;
00421     pi2 = 2.0*pi;
00422     pi2Ts = pi2*Ts;
00423     piDiv2 = pi/2.0;
00424     rad2deg = 180.0/pi;
00425 }
00426 
00427 void GPA::assignFilterStorage()
00428 {
00429     sU = (double*)malloc(3*sizeof(double));
00430     sY = (double*)malloc(3*sizeof(double));   
00431 }
00432 
00433 void GPA::fexcDesLogspace(double fMin, double fMax, int NfexcDes)
00434 {
00435     // calculate logarithmic spaced frequency points
00436     double Gain = log10(fMax/fMin)/((double)NfexcDes - 1.0);
00437     double expon = 0.0;
00438     for(int i = 0; i < NfexcDes; i++) {
00439         fexcDes[i] = fMin*pow(10.0, expon);
00440         expon += Gain;
00441     }
00442 }
00443 
00444 void GPA::calcGPAmeasPara(double fexcDes_i)
00445 {
00446     // Nmeas has to be an integer
00447     Nper = NperMin;
00448     Nmeas = (int)floor((double)Nper/fexcDes_i/Ts + 0.5);
00449     //  secure that the minimal number of measurements is fullfilled
00450     int Ndelta = NmeasMin - Nmeas;
00451     if(Ndelta > 0) {
00452         Nper = (int)ceil((double)NmeasMin*fexcDes_i*Ts);
00453         Nmeas = (int)floor((double)Nper/fexcDes_i/Ts + 0.5);
00454     }
00455     // evaluating reachable frequency
00456     fexc = (double)Nper/(double)Nmeas/Ts;
00457 }
00458 
00459 void GPA::calcGPAsweepPara()
00460 {
00461     // calculate linear frequency sweep parameters
00462     double ksta = ceil(Ts*(double)Nsweep/2.0*(fexc + fexcPast));
00463     Nsweep = (int)floor(2.0*ksta/Ts/(fexc + fexcPast) + 0.5);
00464     bfexc = 2.0*ksta/Ts/(double)Nsweep - fexc;
00465     afexc = (fexc - bfexc)/((double)Nsweep + 1.0);
00466     aAexc = (Aexc - AexcPast)/((double)Nsweep + 1.0);
00467     bAexc = AexcPast;
00468 }
00469 
00470 double GPA::wrapAngle(double angle)
00471 {
00472     // wrap angle from (-2pi,2pi) into (-pi,pi)
00473     if(abs(angle) > pi) angle -= copysign(-pi2, angle); // -1*sign(angle)*2*pi + angle;
00474     return angle;
00475 }
00476 
00477 void GPA::printLongLine()
00478 {
00479     printf("-------------------------------------------------------------------------------------------------------\r\n");
00480 }
00481 
00482 // -----------------------------------------------------------------------------
00483 //      public functions
00484 // -----------------------------------------------------------------------------     
00485 
00486 void GPA::printGPAfexcDes()
00487 {
00488     printLine();
00489     for(int i = 0; i < NfexcDes; i++) {
00490         printf("%9.4f\r\n", (float)fexcDes[i]);
00491     }
00492 }
00493 
00494 void GPA::printGPAmeasPara()
00495 {
00496     printLine();
00497     printf(" fexcDes[Hz]   fexc[Hz]     Aexc      Nmeas   Nper  Nsweep\r\n");
00498     printLine();
00499     for(int i = 0; i < NfexcDes; i++) {
00500         calcGPAmeasPara(fexcDes[i]);
00501         if(fexc == fexcPast || fexc >= fnyq) {
00502             fexc = 0.0;
00503             Aexc = 0.0;
00504             Nmeas = 0;
00505             Nper = 0;
00506             Nsweep = 0;
00507             afexc = 0.0;
00508             bfexc = 0.0;
00509             aAexc = 0.0;
00510             bAexc = 0.0;
00511             
00512         } else {
00513             Aexc = aAexcDes/fexc + bAexcDes;
00514             if(Nsweep > 0)  calcGPAsweepPara();
00515             else{
00516                 afexc = 0.0;
00517                 bfexc = 0.0;
00518                 aAexc = 0.0;
00519                 bAexc = 0.0; 
00520             }
00521             fexcPast = fexc;
00522             AexcPast = Aexc;
00523         }
00524         NmeasTotal += Nmeas;
00525         NmeasTotal += Nsweep;
00526         printf("%11.4e %12.4e %10.3e %7i %6i %7i\r\n", (float)fexcDes[i], (float)fexc, (float)Aexc, Nmeas, Nper, Nsweep);
00527         Nsweep = NsweepMin;
00528     }
00529     printGPAmeasTime();
00530     reset();
00531 }
00532 
00533 void GPA::printFullGPAmeasPara()
00534 {
00535     printLongLine();
00536     printf(" fexcDes[Hz]   fexc[Hz]     Aexc      Nmeas   Nper  Nsweep    afexc      bfexc      aAexc      bAexc\r\n");
00537     printLongLine();
00538     for(int i = 0; i < NfexcDes; i++) {
00539         calcGPAmeasPara(fexcDes[i]);
00540         if(fexc == fexcPast || fexc >= fnyq) {
00541             fexc = 0.0;
00542             Aexc = 0.0;
00543             Nmeas = 0;
00544             Nper = 0;
00545             Nsweep = 0;
00546             afexc = 0.0;
00547             bfexc = 0.0;
00548             aAexc = 0.0;
00549             bAexc = 0.0;
00550             
00551         } else {
00552             Aexc = aAexcDes/fexc + bAexcDes;
00553             if(Nsweep > 0)  calcGPAsweepPara();
00554             else{
00555                 afexc = 0.0;
00556                 bfexc = 0.0;
00557                 aAexc = 0.0;
00558                 bAexc = 0.0; 
00559             }
00560             fexcPast = fexc;
00561             AexcPast = Aexc;
00562         }
00563         NmeasTotal += Nmeas;
00564         NmeasTotal += Nsweep;
00565         printf("%11.4e %12.4e %10.3e %7i %6i %7i %10.3e %10.3e %10.3e %10.3e\r\n", (float)fexcDes[i], (float)fexc, (float)Aexc, Nmeas, Nper, Nsweep, (float)afexc, (float)bfexc, (float)aAexc, (float)bAexc);
00566         Nsweep = NsweepMin;
00567     }
00568     printGPAmeasTime();
00569     reset();
00570 }
00571 
00572 void GPA::printGPAmeasTime()
00573 {
00574     printLine();
00575     printf(" Number of data points :  %11i\r\n", NmeasTotal);
00576     printf(" Measurment time in sec: %12.2f\r\n", (float)((double)NmeasTotal*Ts));
00577 }
00578 
00579 void GPA::printNfexcDes()
00580 {
00581     printLine();
00582     printf("  Number of frequancy points:  %3i\r\n", NfexcDes);
00583 }
00584 
00585 void GPA::printLine()
00586 {
00587     printf("--------------------------------------------------------------------------------\r\n");
00588 }
00589 
00590 GPA::gpadata_t GPA::getGPAdata()
00591 {
00592     return gpaData;
00593 }