/* PZT_poze_graph.c implements various POZE drawing functions. Using the minimal graphic toolkit 'gt' which can set and get RGB-pixels, draw lines, etc. and may give output in compressed JPEG format. PoZeTools v 0.53, december 19, 2006. Copyright (c) 2002-2006 Pieter Suurmond. Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. Any person wishing to distribute modifications to the Software is requested to send the modifications to the original developer so that they can be incorporated into the canonical version. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND ON INFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include #include #include #include #include "gt.h" /* Minimal JPEG graphics toolkit. */ #include "PZT_complex_math.h" /* Complex types and arithmetic. */ #include "PZT_complex_text.h" /* For writePolar() and cutTrailingZeroesAndDot(). */ #include "PZT_poze.h" /* POZE object containing frequency domain data. */ #include "PZT_poze_graph.h" #include "PZT_graph.h" /* For findScaleMarks(). */ #define kPOZE_COLWIDTH_PIX (200) /* Used in dump_GZP() and POZE_draw_z_plane(). */ /*-----------------------------------------------------------------------------*/ static void dump_GZP(RGBPICp pic, int x_left, /* 2 absolute coordinates. */ int y_top, int x_right, int y_bot, const POZE* pz, char* buf) /* Alloc enough elsewhere. */ { int n, smart, x, y; complexPolar pp; char *str, *plur; (void)x_right; /* Surpress 'unused parameter' warning. */ x = x_left; y = y_top; str = buf + sprintf(buf, "G = "); writePolar(&pz->gain, 4, str); /* Gain already in polar format. */ gtDrawString(pic, x, y, 0,0,0, 1, buf); /* Choose 4 decimals precision. */ y += 16; if (pz->zero.number) { if (pz->zero.number > 1) plur = "s:"; else plur = ":"; } else plur = "."; sprintf(buf, "%d Zero%s", pz->zero.number, plur); gtDrawString(pic, x, y, 0,0,0, 1, buf); y += 2; x += 12; /* Indent 2 spaces. */ for (n=0; n < pz->zero.number; n++) /* Print all zeros: */ { y += 12; if (y >= y_bot) { x += kPOZE_COLWIDTH_PIX; y = y_top; } cartesianToPolar(&pz->zero.cartesian[n], &pp); str = buf + sprintf(buf, "Z%d = ", n); writePolar(&pp, 4, str); /* 4 decimals. */ gtDrawString(pic, x, y, 0,0,0, 1, buf); /* 1=over, 0=transparant. */ } y += 16; x -= 12; /* Back 2 spaces. */ if (pz->pole.number) smart = 18; /* Alleen alvast tabuleren als er nog komt.*/ else smart = 0; if (y >= (y_bot - smart)) { x += kPOZE_COLWIDTH_PIX; y = y_top; } if (pz->pole.number) { if (pz->pole.number > 1) plur = "s:"; else plur = ":"; } else plur = "."; sprintf(buf, "%d Pole%s", pz->pole.number, plur); gtDrawString(pic, x, y, 0,0,0, 1, buf); y += 2; x += 12; /* Indent 2 spaces. */ for (n=0; n < pz->pole.number; n++) /* Print all poles: */ { y += 12; if (y >= y_bot) { x += kPOZE_COLWIDTH_PIX; y = y_top; } cartesianToPolar(&pz->pole.cartesian[n], &pp); str = buf + sprintf(buf, "P%d = ", n); writePolar(&pp, 4, str); /* 4 decimals. */ gtDrawString(pic, x, y, 0,0,0, 1, buf); /* 1=over, 0=transparant. */ } } /*---------------------------------------------------------------------------*/ static int POZE_draw_z_plane(RGBPICp pic, /* First 2 args not NULL! */ const POZE* pz, int xstart, int ystart, int width, int height) { int n, m, x, y, x2, y2, smallest, autoDepth, rr, xcenter, ycenter; double fact, angle, unitPixels, xx, yy, larger, hypq, biggestHypq = 0.5; /* Max auto-zoom. */ long occupied; char buffy[128]; unsigned char Rm, Gm, Bm, red; if (!pz) return 1; if (height < 150) return 2; /* columns needed... */ rr = kPOZE_COLWIDTH_PIX * (1 + ((pz->zero.number + pz->pole.number + 3) * 12 / (height-2))); width -= rr; if (width < 150) return 3; /* (3 lines extra, font-height=12) */ /*------------------------------------ Print values G P Z. --------*/ dump_GZP(pic, xstart + width + 1, /* Top left. */ ystart + 1, xstart + width + rr - 1, /* Bottom right. */ ystart + height - 1, pz, buffy); /*-------------------------------------------------------------*/ xcenter = ((1 + width) >> 1) + xstart, /* Calculate */ ycenter = ((1 + height) >> 1) + ystart; /* relative Origin */ /*------------------------------------ BACKROUND ----------------------*/ x = xstart + width - 6; x2 = xstart + 3; /* Gray-blue unit circle and axes. */ gtDrawLine(pic, x, ycenter, x2, ycenter, 160,160,255); /* Real axis. */ gtDrawString(pic, x-24, ycenter+1, 96,96,255, 0, "Real"); /* 1=over, 0=transparant. */ y = ystart + 1; y2 = ystart + height - 2; gtDrawLine(pic, xcenter, y, xcenter, y2, 160,160,255); /* Imaginary axis. */ gtDrawString(pic, xcenter+2, y, 96,96,255, 0, "Imaginary"); /* Big circle. */ if (width < height) smallest = width; /* Take smallest. */ else smallest = height; /*--------------------------------------- Autozoom: --------------*/ for (n=0; n < pz->zero.number; n++) /* Find largest magnitude. */ { hypq = pz->zero.cartesian[n].re * pz->zero.cartesian[n].re + pz->zero.cartesian[n].im * pz->zero.cartesian[n].im; if (hypq > biggestHypq) biggestHypq = hypq; } for (n=0; n < pz->pole.number; n++) { hypq = pz->pole.cartesian[n].re * pz->pole.cartesian[n].re + pz->pole.cartesian[n].im * pz->pole.cartesian[n].im; if (hypq > biggestHypq) biggestHypq = hypq; } unitPixels = 0.42 * (double)smallest / sqrt(biggestHypq); autoDepth = (int)(7.0 * unitPixels); /* Take 84% of smallest dim. */ fact = kPZT_twoPi / (double)autoDepth; for (n=0; n= 0.0) xx += 0.5; else xx -= 0.5; if (yy >= 0.0) yy += 0.5; else yy -= 0.5; x = xcenter + (int)xx; y = ycenter + (int)yy; if ((y > ystart) && (y < ystart + height) && (x > xstart) && (x < xstart + width)) gtSetPixel(pic, x, y, 160,160,255); } /*------------------------------------------------------------*/ for (n=0; n < pz->zero.number; n++) /* Draw all zeros: */ { xx = unitPixels * pz->zero.cartesian[n].re; yy = unitPixels * pz->zero.cartesian[n].im; if (xx >= 0.0) xx += 0.5; else xx -= 0.5; if (yy >= 0.0) yy += 0.5; else yy -= 0.5; x = xcenter + (int)xx; y = ycenter - (int)yy; larger = 2.8; /* Find out how many pixels are occupied */ do { /* by black. If many, grow the circle... */ larger += 4.0; occupied = 0L; autoDepth = (int)(5.71 * larger); fact = kPZT_twoPi / (double)autoDepth; for (m=0; m>= 5 /* was 7 !!!!!*/; /* half unit circle. */ } while ((occupied >= (long)autoDepth) && (larger <= 0.5 * unitPixels)); if (larger > 0.5 * unitPixels) red = 224; /* Mark as red when we didn't grow. */ else red = 0; for (m=0; mpole.number; n++) /* Draw all poles: */ { x = xcenter + (int)(unitPixels * pz->pole.cartesian[n].re); y = ycenter - (int)(unitPixels * pz->pole.cartesian[n].im); x2 = x - 5; y2 = y - 5; /* No GROW-ROTATE ALGORITHM YET. */ x += 5; y += 5; gtDrawLine(pic, x, y, x2, y2, 0,0,0); /* X */ gtDrawLine(pic, x, y-10, x2, y2+10, 0,0,0); } return 0; } /*-----------------------------------------------------------------------------*/ int POZE_jpeg_z_plane (const POZE *pz, int width, int height, const char *title, int q, /* JPEG quality. */ FILE* fp) /* Destination. */ { RGBPICp pic; int e, hd; char buffy[128]; /* BACKROUNDS, AXES, ETC: */ if (gtCreateMemRGB(&pic, width, height, 255,255,255)) /* White background. */ return 1; /* 1 means cannot create. */ if (title) { gtDrawString(pic, 2,0, 128,128,128, /* RGB colour. */ 1, title); /* 1=over, 0=transparant. */ hd = 10; /* Top margin reserved for title. */ } else hd = 0; /* startX: startY: width: height: */ e = POZE_draw_z_plane(pic, pz, 0, hd, width, height - hd); if (e) { sprintf(buffy, "POZE_draw_z_plane() = %d!", e); gtDrawString(pic, 3, 3, 255,0,0, 1, buffy); e = 2; /* 2 means error during drawing of poles and zeros. */ } else { e = gtOutputMemRGB(pic, fp, q); if (e) e = 3; /* 3 means could not save. */ } gtDestroyMemRGB(&pic); return e; } /*-------------------------------------------------------------*/ /* Draws a big PI. Specify topleft position. */ static void draw_large_pi(RGBPICp pic, int x, int y) { /* \ ---- \ */ gtDrawLine(pic, x, y, x+2, y+2, 0,0,0); gtDrawLine(pic, x+14, y+2, x+2, y+2, 0,0,0); gtDrawLine(pic, x+14, y+2, x+16, y+4, 0,0,0); /* || */ gtDrawLine(pic, x+5, y+2, x+5, y+17, 0,0,0); gtDrawLine(pic, x+11, y+2, x+11, y+17, 0,0,0); } /*-----------------------------------------------------------------------------*/ /* Draws the definition (formula) of the transfer function H(w) from the complex Z-plane. There is another routine for the definition from the time domain. Specify upper left corner within the picture. It always takes a fixed width and height: 140 x 100. */ static int draw_H_def_f(RGBPICp pic, /* First 2 args not NULL! */ const POZE *pz, int xstart, int ystart) { int n, x, y; char buffy[16]; if (!pic) return 1; if (!pz) return 2; x = xstart + 1; y = ystart + 48; /* 0,0,0 means select black. */ gtDrawString(pic, x,y, 0,0,0, 1, "H"); /* 1=over, 0=transparant. */ x += 6; y += 2; gtDrawString(pic, x,y, 0,0,0, 1, "(w)"); x += 20; y -= 2; gtDrawString(pic, x,y, 0,0,0, 1, "= G"); x += 22; y += 5; /* Divider-line. */ gtDrawLine(pic, x, y, x + 84, y, 0,0,0); x += 12; y -= 12; gtDrawString(pic, x, y, 0,0,0, 1, "m=0"); y -= 29; n = 3 * (sprintf(buffy, "m<%d", pz->zero.number) - 3); x -= n; /* center. */ gtDrawString(pic, x, y, 0,0,0, 1, buffy); x += n; y += 11; draw_large_pi(pic, x, y); x += 21; y += 6; gtDrawString(pic, x, y, 0,0,0, 1, "e"); x += 6; y -= 6; gtDrawString(pic, x, y, 0,0,0, 1, "jw"); x += 14; y += 6; gtDrawString(pic, x, y, 0,0,0, 1, "- Z"); x += 19; y += 4; gtDrawString(pic, x, y, 0,0,0, 0, "m"); x -= 60; /*--------------------- */ y += 50; gtDrawString(pic, x, y, 0,0,0, 1, "n=0"); y -= 29; n = 3 * (sprintf(buffy, "n<%d", pz->pole.number) - 3); x -= n; /* center. */ gtDrawString(pic, x, y, 0,0,0, 1, buffy); x += n; y += 11; draw_large_pi(pic, x, y); x += 21; y += 6; gtDrawString(pic, x, y, 0,0,0, 1, "e"); x += 6; y -= 6; gtDrawString(pic, x, y, 0,0,0, 1, "jw"); x += 14; y += 6; gtDrawString(pic, x, y, 0,0,0, 1, "- P"); x += 18; y += 4; gtDrawString(pic, x, y, 0,0,0, 0, "n"); return 0; } /*----------------------------------------------------------------------------*/ /* When log_freq_scale is non-zero, always 10 octaves are shown. In case it is zero, full linear frequency range from 0 to pi is shown. There may be cases that we cannot plot correctly: it all has to do with phase unwrapping: when poles or zeros have very near angles, we must be sure that the phase-unwrapper calculates enough points between them. Also the first angle, in case of logarithmic frequency scale, may skip a first pole or zero, so we're unable to correctly unfold phase. v0.53: Added a clipper for time. */ static int POZE_draw_response_freq(RGBPICp pic, /* Trusted. */ const POZE* pz, /* Checked. */ int xstart, /* Absolute */ int ystart, /* offset. */ int width, /* Subsize. */ int height, short log_freq_scale) { static const int y_dim_top = 16, y_text = 14, /* Space between graphs.*/ xbefore = 50, /* left margin. */ xafter = 8, /* right margin. */ numOctaves = 10; /* Constant for now. */ int y_lin_top, y_lin_bottom, y_log_top, y_log_bottom, y_ang_top, y_ang_bottom, y_tim_top, y_tim_bottom, eff_h, x_begin, x_end, x_pix, n, m, x, y, div, RG; double h, pixelDensity, startFreq, fact, lconst, xx, yy, step, steplin, botMagLin, topMagLin, topMagLog, botMagLog, botAng, topAng, botTim, topTim, unitYPixels_lin, unitYPixels_log, unitYPixels_ang, unitYPixels_tim, angle, prev_angle, prev_omega; /* Used for phase unwrapping. */ char buffy[128]; complexPolar z_p, H_p; complexCartesian z_c; response_range r_range; if (!pz) return 1; eff_h = (height - (y_text * 8)) >> 2; /* Four plots. */ if (eff_h < 32) return 2;/* Too less pixels left. */ y_lin_top = ystart + y_text + 1; y_lin_bottom = y_lin_top + eff_h; y_log_top = y_lin_bottom + (y_text << 1) + 1; y_log_bottom = y_log_top + eff_h; y_ang_top = y_log_bottom + (y_text << 1); y_ang_bottom = y_ang_top + eff_h; y_tim_top = y_ang_bottom + (y_text << 1); y_tim_bottom = y_tim_top + eff_h; /*------------------------------------------ Amplitude: -------------------*/ if (log_freq_scale) startFreq = kPZT_onePi / pow(2.0, (double)numOctaves); else startFreq = 0.0; /* startFreq also used to init z_p.ang, far below. */ POZE_response_range(pz, &r_range); /* ALWAYS USES 0 AS STARTFREQ! */ if ((r_range.max_magnitude <= 0.0001) || (r_range.max_magnitude >= 10000.0)) return 3; /* That's too small / large. */ /*------------------------------------ CHOOSE HORIZONTAL DIVISIONS: --------*/ x_begin = xstart + xbefore; x_end = xstart + width - xafter; x_pix = x_end - x_begin; if (x_pix < 100) return 4; /* Too less pixels left. */ if (log_freq_scale) { if (x_pix < 200) div = 2; /* Keep in mind we draw 10 octaves. */ else if (x_pix < 500) div = 5; else if (x_pix < 1000) div = 10; else div = 20; } else { if (x_pix < 200) div = 4; /* Linear frequency scale. */ else if (x_pix < 300) div = 5; else if (x_pix < 450) div = 8; else if (x_pix < 600) div = 10; else if (x_pix < 1200) div = 20; else div = 40; } /*------------------------------------ DRAW VERTICAL LINES: -------------*/ for (n = 0; n <= div; n++) { if (n) RG = 196; /* Somewhat lighter. */ else RG = 160; fact = (double)n / (double)div; x = x_begin + (int)(0.5 + ((double)x_pix * fact)); gtDrawLine(pic, x, y_lin_top-1, x, y_lin_bottom+1, RG,RG,255); /* Vertical lin magnitude axis. */ gtDrawLine(pic, x, y_log_top-1, x, y_log_bottom+1, RG,RG,255); /* Vertical log magnitude axis. */ gtDrawLine(pic, x, y_ang_top-1, x, y_ang_bottom+1, RG,RG,255); /* Vertical phase axis. */ gtDrawLine(pic, x, y_tim_top-1, x, y_tim_bottom+1, RG,RG,255); /* Vertical time axis. */ if (log_freq_scale) /* Always 10 octaves: */ fact = pow(2.0, (double)(numOctaves*(n-div)) / (double)div); sprintf(buffy, "%.4f", fact); m = cutTrailingZeroesAndDot(buffy); /* (m=strlen.) */ if (strcmp(buffy, "0")) /* Append pi char when not zero. */ { if (strcmp(buffy, "1")) sprintf(buffy + m, "~"); else /* One pi needs no 1, replace. */ sprintf(buffy, "~"); } x -= (3 * strlen(buffy)) - 1; gtDrawString(pic, x, y_lin_bottom+2, 96,96,255, 0, buffy); gtDrawString(pic, x, y_log_bottom+2, 96,96,255, 0, buffy); gtDrawString(pic, x, y_ang_bottom+2, 96,96,255, 0, buffy); gtDrawString(pic, x, y_tim_bottom+2, 96,96,255, 0, buffy); } /*-------------------------------- linear amplitude ---------------------------------*/ gtDrawString(pic, xstart + 4, y_lin_top - y_dim_top, 0,0,0, 0, "|H |"); gtDrawString(pic, xstart + 4, y_lin_top - y_dim_top + 1, 0,0,0, 0, " (w)"); arrow(pic, xstart + 10, y_lin_top - 2, 1); arrow(pic, (xstart + x_end) >> 1, y_lin_top - 8, 0); gtDrawString(pic, ((xstart + x_end) >> 1) + 4, y_lin_top - 14, 0,0,0, 0, "w"); botMagLin = r_range.min_magnitude; topMagLin = r_range.max_magnitude; n = findScaleMarks (eff_h, &botMagLin, &topMagLin, &step); if (n) return 1000 + n; unitYPixels_lin = (double)eff_h / (topMagLin - botMagLin); for (yy = botMagLin; yy < (topMagLin + kPZT_epsilon); yy += step) { if ((yy != botMagLin) && (yy != topMagLin)) RG = 196; /* Somewhat lighter. */ else RG = 160; y = y_lin_bottom - (int)(0.5 + (unitYPixels_lin * (yy - botMagLin))); gtDrawLine(pic, x_begin - 1, y, x_end + 1, y, RG,RG,255); /* Horiz. freq. axis. */ sprintf(buffy, "%.6f", yy); /* (m=strlen.) */ m = cutTrailingZeroesAndDot(buffy); y -= 6; if (fabs(yy - topMagLin) < kPZT_epsilon) y++; gtDrawString(pic, x_begin - 2 - (6 * m), y, 96,96,255, 0, buffy); /* 1=over, 0=transparant. */ } /*--------------------------- log amplitude ---------------------------------*/ gtDrawString(pic, xstart + 4, y_log_bottom-18, 96,96,255, 0, "dB"); gtDrawString(pic, xstart + 4, y_log_top - y_dim_top - 1, 0,0,0, 0, " 10"); gtDrawString(pic, xstart + 4, y_log_top - y_dim_top, 0,0,0, 0, "20log"); gtDrawString(pic, xstart + 4, y_log_top - y_dim_top + 1, 0,0,0, 0, " (|H |)"); gtDrawString(pic, xstart + 4, y_log_top - y_dim_top + 2, 0,0,0, 0, " (w)"); arrow(pic, xstart + 10, y_log_top - 2, 1); arrow(pic, (xstart + x_end) >> 1, y_log_top - 8, 0); gtDrawString(pic, ((xstart + x_end) >> 1) + 4, y_log_top - 14, 0,0,0, 0, "w"); topMagLog = 20.0 * log10(r_range.max_magnitude); /* Dynamic range <= 120 dB. */ if (r_range.min_magnitude <= (r_range.max_magnitude / 1000000.0)) botMagLog = topMagLog - 96.0; /* 120.0 dB corresponds to a ratio of 1 million. */ else botMagLog = 20.0 * log10(r_range.min_magnitude); n = findScaleMarks (eff_h, &botMagLog, &topMagLog, &step); if (n) return 2000 + n; unitYPixels_log = (double)eff_h / (topMagLog - botMagLog); for (yy = botMagLog; yy < (topMagLog + kPZT_epsilon); yy += step) { if (fabs(yy) < kPZT_epsilon) RG = 160; /* 0 dB somewhat darker. */ else RG = 196; y = y_log_bottom - (int)(0.5 + (unitYPixels_log * (yy - botMagLog))); gtDrawLine(pic, x_begin-1, y, x_end+1, y, RG,RG,255); /* Horizontal freq. axis. */ sprintf(buffy, "%.6f", yy); m = cutTrailingZeroesAndDot(buffy); y -= 6; if (fabs(yy - topMagLog) < kPZT_epsilon) y++; gtDrawString(pic, x_begin-2-(6*m), y, 96,96,255, 0, buffy); /* 1=over, 0=transparant. */ } /*--------------------------- ANGLE ------------------------------------------*/ gtDrawString(pic, xstart + 4, y_ang_top - y_dim_top, 0,0,0, 0, "> 1, y_ang_top - 8, 0); gtDrawString(pic, ((xstart + x_end) >> 1) + 4, y_ang_top - 14, 0,0,0, 0, "w"); botAng = r_range.min_angle / kPZT_onePi; /* Later, at plotting, H_p.ang */ topAng = r_range.max_angle / kPZT_onePi; /* is also divided by PI. */ n = findScaleMarks (eff_h, &botAng, &topAng, &step); if (n) return 3000 + n; unitYPixels_ang = (double)eff_h / (topAng - botAng); for (yy = botAng; yy < (topAng + kPZT_epsilon); yy += step) { h = fabs(yy / kPZT_twoPi); if (fabs(h - trunc(0.5 + h)) < kPZT_epsilon) RG = 160; else RG = 196; x = x_begin - 1; /* left. */ y = y_ang_bottom - (int)(0.5 + (unitYPixels_ang * (yy - botAng))); gtDrawLine(pic, x, y, x_end+1, y, RG,RG,255); /* Horizontal freq. axis. */ sprintf(buffy, "%.4f", yy); m = cutTrailingZeroesAndDot(buffy); if (strcmp(buffy, "0")) /* Append pi char when not zero. */ { if (!strcmp(buffy, "1")) /* Plus or minus pi needs no 1, replace. */ sprintf(buffy, "~"); else if (!strcmp(buffy, "-1")) /* Plus or minus pi needs no 1, replace. */ sprintf(buffy, "-~"); else sprintf(buffy + m, "~"); } x -= 2 + (6 * strlen(buffy)); y -= 6; if (fabs(yy - topAng) < kPZT_epsilon) y++; gtDrawString(pic, x, y, 96,96,255, 0, buffy); /* Somewhat darker again. */ } /*--------------------------- TIME - DELAY (angle/-omega): --------------------------*/ gtDrawString(pic, xstart + 4, y_tim_bottom - 18, 96,96,255, 0, "Ts"); gtDrawString(pic, xstart + 4, y_tim_top - y_dim_top, 0,0,0, 0, "-d> 1, y_tim_top - 8, 0); gtDrawString(pic, ((xstart + x_end) >> 1) + 4, y_tim_top - 14, 0,0,0, 0, "w"); botTim = r_range.min_time; topTim = r_range.max_time; n = findScaleMarks (eff_h, &botTim, &topTim, &step); if (n) return 4000 + n; unitYPixels_tim = (double)eff_h / (topTim - botTim); for (yy = botTim; yy < (topTim + kPZT_epsilon); yy += step) { if (fabs(yy) < kPZT_epsilon) RG = 160; else RG = 196; y = y_tim_bottom - (int)(0.5 + (unitYPixels_tim * (yy - botTim))); x = x_begin - 1; /* left. */ gtDrawLine(pic, x, y, x_end+1, y, RG,RG,255); /* Horizontal freq. axis. */ sprintf(buffy, "%.4f", yy); m = cutTrailingZeroesAndDot(buffy); x -= 2 + (6 * m); y -= 6; if (fabs(yy - topTim) < kPZT_epsilon) y++; gtDrawString(pic, x, y, 96,96,255, 0, buffy); } /*-------------------------------- THE ACTUAL PLOTTING -------------------------------*/ pixelDensity = 2.0 * (double)x_pix * (double)eff_h; /* Adapt to size. */ steplin = kPZT_onePi / pixelDensity; /* Needed even for log freq scale. */ if (log_freq_scale) { fact = (double)(x_pix) / log(pow(2.0, (double)numOctaves)); lconst = log(pow(2.0,(double)numOctaves) / kPZT_onePi); step = pow(2.0, (double)numOctaves / pixelDensity); } /* Calculated y-pixels per x-pixel. */ else { lconst = 1.0; /* 1.0 = dummy value to satisfy compiler (surpress uninitialised warrning). */ fact = (double)(x_pix) / kPZT_onePi; } prev_angle = r_range.start_angle; /* Calculated by POZE_response_range() far above. */ z_p.ang = prev_omega = r_range.start_omega; z_p.mag = 1.0; x = 0; /* Just to satisfy the compiler. */ while (z_p.ang <= kPZT_onePi) { if (log_freq_scale && (z_p.ang >= startFreq)) z_p.ang *= step; else z_p.ang += steplin; polarToCartesian(&z_p, &z_c); POZE_response(pz, &z_c, &H_p); /* Polar result. */ if (z_p.ang >= startFreq) { if (log_freq_scale) xx = log(z_p.ang) + lconst; /* Log. */ else xx = z_p.ang; /* Lin. */ x = x_begin + (int)(0.5 + (fact * xx)); /*------------- plot lin magnitude: -------------------------------------*/ y = y_lin_bottom - (int)(0.5 + (unitYPixels_lin * (H_p.mag - botMagLin))); gtSetPixel(pic, x, y, 0,0,0); /*------------- plot log magnitude: -------------------------------------*/ if (H_p.mag > kPZT_epsilon) { y = y_log_bottom - (int) (0.5 + (unitYPixels_log * (20.0 * log10(H_p.mag) - botMagLog))); if (y <= y_log_bottom) /* We may exceed 120 dB. */ gtSetPixel(pic, x, y, 0,0,0); } } /*------------- calculate phase: --------------------------------------------*/ if (H_p.mag >= kPZT_epsilon) { angle = phaseUnwrap(prev_angle, H_p.ang); /* Always unwrap. */ if (z_p.ang >= startFreq) { y = y_ang_bottom - (int) (0.5 + (unitYPixels_ang * ((angle / kPZT_onePi) - botAng))); gtSetPixel(pic, x, y, 0,0,0); /*------------- plot time-displacement: -----------------------------*/ y = y_tim_bottom - (int) (0.5 + (unitYPixels_tim * ( ((prev_angle - angle) / (z_p.ang - prev_omega)) - botTim))); /* We need to clip if POZE_response_range() excludes points that are too large... */ if ((y <= y_tim_bottom) && (y >= y_tim_top)) gtSetPixel(pic, x, y, 0,0,0); /* We may skip points around zeros! */ } prev_omega = z_p.ang; prev_angle = angle; } } return 0; } /*------------------------------------------------------------------------*/ int POZE_jpeg_response_freq(const POZE* pz, /* Tested in subfunc. */ int width, int height, const char* title, /* Title may be NULL. */ int q, /* JPEG quality. */ FILE* fp) /* Destination. */ { RGBPICp pic; /* Generates jpeg-image in the current directory. */ int hd; short e; char buffy[128]; if (gtCreateMemRGB(&pic, width, height, 255, 255, 255)) /* White bg. */ return 1; if (title) { gtDrawString(pic, 2, 0, 128,128,128, /* R, G, B (grey). */ 1, title); /* 1=over, 0=transparant. */ hd = 11; /* Top margin reserved for title. */ } else hd = 0; /* Argument 'z_plane' is tested here. */ e = POZE_draw_response_freq(pic, pz, /* Pass P,Z and G. */ 0, hd, /* startx, starty. */ width, height-hd, /* width, height. */ 0); /* scale settings. */ if (e) { sprintf(buffy, "POZE_draw_response_freq() = %d!", e); gtDrawString(pic, width >> 3, height >> 1, 255, 255, 255, 1, buffy); } e = gtOutputMemRGB(pic, fp, q); gtDestroyMemRGB(&pic); return (int)e; } /*-------------------------------------------------------------------*/ int POZE_jpeg_combined_freq(const POZE* pz, int log_freq, int width, int height, const char* title, int q, /* JPEG quality. */ FILE* fp) /* Destination. */ { RGBPICp pic; int n, hd, h1, h2; char buffy[256]; short e, ee; /*-------------------------------- BACKROUNDS, AXES, ETC: ----------------*/ if (gtCreateMemRGB(&pic, width, height, 255, 255, 255)) /* White bg. */ return 1; if (title) { hd = 12; /* Top margin for title. */ gtDrawString(pic, 2, 0, 128,128,128, 1, title); /* 1=over,0=transprt. */ n = 0; if (!POZE_stable(pz)) /* (Results are all 0 or 1.) */ n += sprintf(buffy, "un"); n += sprintf(buffy + n, "stable "); if (!POZE_causal(pz)) n += sprintf(buffy + n, "non-"); n += sprintf(buffy + n, "causal "); if (POZE_normalized((POZE*)pz, 0)) /* Doit = 0 means we */ n += sprintf(buffy + n, "normalized "); /* don't touch *pz. */ if (POZE_fir(pz)) n += sprintf(buffy + n, "FIR"); else n += sprintf(buffy + n, "IIR"); buffy[0] -= 32; /* Upcase first letter. */ gtDrawString(pic, (1 + width - (6 * n)) >> 1, 0, 0,0,0, 1, buffy); } else hd = 0; h1 = (1 + height - hd) / 4; /* Top part 25%. */ n = 12 * (4 + pz->zero.number + pz->pole.number); if (n > h1) h1 = n; h2 = (height - hd) - h1; /* Bottom remaining 75%. */ /*-------------------- DRAW POLES AND ZEROES (and gain): ------------*/ ee = POZE_draw_z_plane(pic, pz, /* Pass P,Z and G, */ 0, hd, /* Startx, starty. */ width-140, h1); /* Width, height. */ if (ee) { sprintf(buffy, "ERROR POZE_draw_z_plane() = %d!", ee); gtDrawString(pic, 16, h1 >> 1, 255,0,0, 1, buffy); } /*--------------------- DRAW FORMULA ON RIGHT-SIDE: -------------------*/ ee = draw_H_def_f(pic, pz, width-140, (h1>>1)-80); if (ee) /* ..... later: center within left Y-space !!! */ { sprintf(buffy, "ERROR draw_H_def_f() = %d!", ee); gtDrawString(pic, width-138, 50, 255,0,0, 1, buffy); } /*--------------------- DRAW TRANSFER GRAPHS: -------------------------*/ ee = POZE_draw_response_freq(pic, pz, /* Pass P,Z and G, */ 0, h1 + hd, /* Startx, starty. */ width, h2, /* Width, height. */ log_freq); /* Scale settings. */ if (ee) { sprintf(buffy, "ERROR POZE_draw_response_freq() = %d!", ee); gtDrawString(pic, 16, h1 + (h2 >> 1), 255,0,0, 1, buffy); } /*--------------------- SAVE PICTURE TO DISK: -------------------------*/ e = gtOutputMemRGB(pic, fp, q); gtDestroyMemRGB(&pic); return (int)e; }