share/z_transform/z_transform.mac

   1 /* z-transform code
   2  * Copyright 2007 by Robert Dodier
   3  * I release this work under terms of GNU General Public License
   4  *
   5  * Summary. Z-transforms for various special cases are implemented as
   6  * pattern-matching rules. Given an expression like z_transform(foo(n, z), n, z),
   7  * for most rules it's necessary to have the variables n and z in hand before
   8  * looking at foo(n, z). However, the Maxima pattern matcher does not have
   9  * backtracking, so one rule (r0 below) captures n and z, then all other
  10  * rules are applied.
  11  *
  12  * The inverse transform isn't implemented yet. For some inspiration, see:
  13  * http://ece.citadel.edu/barsanti/elec407/L3%20Inverse%20Z%20Transforms.pdf
  14  * and: http://www.reduce-algebra.com/docs/ztrans.pdf
  15  */
  16
  17 put ('z_transform, true, 'present);
  18
  19 apply_z_transform (e) := apply1 (e, r0);
  20
  21 matchdeclare ([nn, zz], symbolp);
  22 matchdeclare (aa, all);
  23
  24 defrule (r0,
  25     z_transform (aa, nn, zz),
  26     block ([nn% : nn, zz% : zz],
  27         apply1
  28            (z_transform (aa, nn, zz),
  29             r913_1a, r913_1b, r913_2a, r913_2b, r913_3a, r913_3b,
  30             r913_4, r913_5, r913_6, r913_7, r913_10, r913_12, r914_4,
  31             r914_6, /* r914_9, */ r914_10, r914_11a, r914_11b, r914_12,
  32             r914_13, r914_14, r914_15a, r914_16)));
  33
  34 /* Some specific transforms.
  35  * Table 9.1.3 at: http://mathfaculty.fullerton.edu/MATHEWS/C2003/ZTRANSFORMINTROMOD.HTML
  36  */
  37
  38 /* (1) delta[n] --> 1
  39  * (have to try kron_delta both ways ... sigh)
  40  */
  41
  42 simp : false;
  43
  44 defrule (r913_1a,
  45     z_transform (kron_delta (nn%, aa), nn%, zz%),
  46     zz^(- aa));
  47
  48 defrule (r913_1b,
  49     z_transform (kron_delta (aa, nn%), nn%, zz%),
  50     zz^(- aa));
  51
  52 simp : true;
  53
  54 /* (2) u[n] --> z/(z - 1) */
  55
  56 defrule (r913_2a,
  57     z_transform (1, nn%, zz%),
  58     zz/(zz - 1));
  59
  60 defrule (r913_2b,
  61     z_transform (unit_step (nn%), nn%, zz%),
  62     zz/(zz - 1));
  63
  64 /* (3) b^n --> z/(z - b) */
  65
  66 matchdeclare (bb, freeof (nn%, zz%));
  67
  68 defrule (r913_3a,
  69     z_transform (bb^nn%, nn%, zz%),
  70     zz/(zz - bb));
  71
  72 defrule (r913_3b,
  73     z_transform (1/(bb^nn%), nn%, zz%),
  74     zz/(zz - 1/bb));
  75
  76 /* (4) b^(n - 1) * u[n - 1] --> 1/(z - b) */
  77
  78 defrule (r913_4,
  79     z_transform (bb^(nn% - 1) * unit_step (nn% - 1), nn%, zz%),
  80     1 / (zz - bb));
  81
  82 /* (5) e^(a*n) --> z/(z - e^a) */
  83
  84 matchdeclare (aa, lambda ([e], e # 0 and freeof (nn%, zz%, e)));
  85
  86 defrule (r913_5,
  87     z_transform (exp (aa * nn%), nn%, zz%),
  88     zz / (zz - exp (aa)));
  89
  90 /* (6) n --> z/(z - 1)^2 */
  91
  92 defrule (r913_6,
  93     z_transform (nn%, nn%, zz%),
  94     zz / (zz - 1)^2);
  95
  96 /* (7) n^2 --> z*(z + 1)/(z - 1)^3 */
  97
  98 defrule (r913_7,
  99     z_transform (nn%^2, nn%, zz%),
 100     zz*(zz + 1) / (zz - 1)^3);
 101
 102 /* (8) b^n*n --> b*z/(z - b)^2
 103  * via (6) + frequency scaling
 104  */
 105
 106 /* (9) e^(a*n)*n --> z*e^a/(z - e^a)^2
 107  * via (6) + complex translation
 108  */
 109
 110 /* (10) sin(a*n) --> sin(a)*z/(z^2 - 2*cos(a)*z + 1) */
 111
 112 defrule (r913_10,
 113     z_transform (sin (aa*nn%), nn%, zz%),
 114     sin(aa)*zz / (zz^2 - 2*cos(aa)*zz + 1));
 115
 116 /* (11) b^n*sin(a*n) --> sin(a)*b*z/(z^2 - 2*cos(a)*b*z + b^2)
 117  * via (10) + frequency scaling
 118  */
 119
 120 /* (12) cos(a*n) --> z*(z - cos(a))/(z^2 - 2*cos(a)*z + 1) */
 121
 122 defrule (r913_12,
 123     z_transform (cos (aa*nn%), nn%, zz%),
 124     zz*(zz - cos(aa)) / (zz^2 - 2*cos(aa)*zz + 1));
 125
 126 /* (13) b^n*cos(a*n) --> z*(z - b*cos(a))/(z^2 - 2*cos(a)*b*z + b^2)
 127  * via (11) + frequency scaling
 128  */
 129
 130
 131 /* General properties.
 132  * Table 9.1.4 at: http://mathfaculty.fullerton.edu/MATHEWS/C2003/ZTRANSFORMINTROMOD.HTML
 133  */
 134
 135 /* (4) u[n - m] --> z^(1 - m)/(z - 1)
 136  * (delayed unit step)
 137  */
 138
 139 matchdeclare (mm, integerp);
 140
 141 defrule (r914_4,
 142     z_transform (unit_step (nn% - mm), nn%, zz%),
 143     zz^(1 - mm) / (zz - 1));
 144
 145 /* (5) x[n - 1]*u[n - 1] --> (1/z)*X(z)
 146  * via (6) w/ m = 1
 147  */
 148
 149 /* (6) x[n - m]*u[n - m] --> z^(-m)*X(z)
 150  * (time delayed shift)
 151  */
 152
 153 defrule (r914_6,
 154     z_transform (aa * unit_step (nn% - mm), nn%, zz%),
 155     z^(-m) * z_transform (subst (nn + mm, nn, aa), nn, zz));
 156
 157 /* (7) x[n + 1] --> z*(X(z) - x[0])
 158  * (8) x[n + 2] --> z^2*(X(z) - x[0] - x[1]*z^(-1))
 159  * via (9) w/ m = 1 and m = 2 respectively
 160  */
 161
 162 /* (9) x[n + m] --> z^m*(X(z) - sum(x[i]*z^(-i), i, 0, m - 1))
 163  * (time forward)
 164  */
 165
 166 /* HMM, NOT SURE HOW TO DO (9) ... FOLLOWING STUFF IS BROKEN */
 167
 168 matchdeclare (mm, lambda ([e], integerp(e) and e > 0));
 169 defmatch (n_plus_m, nn% + mm);
 170
 171 matchdeclare (xxnpm, lambda ([e], n_plus_m (e) # false));
 172
 173 defrule (r914_9,
 174     z_transform (ss, nn%, zz%),
 175     zz^mm * z_transform (subst (nn - mm, nn, ss), nn, zz));
 176
 177 /* (10) e^(a*n)*x[n] --> X(z*e^(-a))
 178  * (complex translation)
 179  */
 180
 181 matchdeclare (nz, lambda ([e], e # 0 and freeof (nn%, zz%, e)));
 182
 183 defrule (r914_10,
 184     z_transform (exp (nz * nn%) * bb, nn%, zz%),
 185     'subst (zz/exp(nz), zz, z_transform (bb, nn, zz)));
 186
 187 /* (11) b^n*x[n] --> X(z/b)
 188  * (frequency scaling)
 189  */
 190
 191 matchdeclare (xx, all);
 192 matchdeclare (bb, freeof (nn%, zz%));
 193
 194 defrule (r914_11a,
 195     z_transform (bb^nn% * xx, nn%, zz%),
 196     'subst (zz/bb, zz, z_transform (xx, nn, zz)));
 197
 198 defrule (r914_11b,
 199     z_transform (1/(bb^nn%) * xx, nn%, zz%),
 200     'subst (zz*bb, zz, z_transform (xx, nn, zz)));
 201
 202 /* (12) n*x[n] --> -z X'(z)
 203  * (differentiation)
 204  */
 205
 206 matchdeclare (aa, all);
 207 matchdeclare (kk, lambda ([e], integerp(e) and e > 0));
 208
 209 defrule (r914_12,
 210     z_transform (aa*nn%^kk, nn%, zz%),
 211     block ([ee : - zz * 'diff (z_transform (aa, nn, zz), zz)],
 212         for i:2 thru kk do ee : - zz * 'diff (ee, zz), ee));
 213
 214 /* (13) (1/n)*x[n] --> - \int X(z)/z dz
 215  * (integration)
 216  */
 217
 218 matchdeclare (uu, lambda ([e], not atom(e) and op(e) = "/" and member (nn%, second(e))));
 219
 220 defrule (r914_13,
 221     z_transform (uu, nn%, zz%),
 222     'integrate (zz^-1 * z_transform (aa, nn, zz), zz));
 223
 224 /* (14) 1/(n + m)*x[n] --> - z^(-m) * \int X(z)/z^(m + 1) dz
 225  * (integration shift)
 226  */
 227
 228 matchdeclare (mm, lambda ([e], integerp(e) and e > 0));
 229
 230 defrule (r914_14,
 231     z_transform (aa/(nn% + mm), nn%, zz%),
 232     - zz^(-mm) * 'integrate (z^-(mm + 1) * z_transform (aa, nn, zz), zz));
 233
 234 /* (15) x[n] (star) y[n] = \sum_{i=0}^n x[i]*y[n - i] --> X(z)*Y(z)
 235  * (discrete time convolution)
 236  */
 237
 238 matchdeclare (cc, lambda ([e], not atom(e) and op(e) = 'convolution));
 239
 240 defrule (r914_15a,
 241     z_transform (cc, nn%, zz%),
 242     block ([a : args(cc)], product (z_transform (a[i], nn, zz), i, 1 , length(a))));
 243
 244 /* ANOTHER RULE R914_15B FOR EXPLICIT 'SUM EXPRESSION WOULD BE NICE ... */
 245
 246 /* (16) \sum_{i=0}^n x[i] --> z/(z - 1)*X(z)
 247  * (convolution with y[n] = 1)
 248  */
 249
 250 matchdeclare (ii, symbolp);
 251
 252 simp : false;
 253
 254 defrule (r914_16,
 255     z_transform ('sum (aa, ii, 0, nn%), nn%, zz%),
 256     zz/(zz - 1) * z_transform (subst (nn, ii, aa), nn, zz));
 257
 258 simp : true;
 259
 260 /* ((17) & (18) -- not transform pairs) */
 261
 262 /* (put linearity declaration last, otherwise messes up rules) */
 263
 264 /* (1) addition
 265  * (2) constant multiple
 266  * (3) linearity
 267  */
 268
 269 declare (z_transform, linear);
 270