zpk2tf and freqz functions implementation (#2988)

* Add Zpk2tf function

* Add zpk2tf function documentation

* Fix

* Add freqz function implementation

* Fix format

* Add freqz function documentation

* Fix

* Fix format

* Update index.d.ts

* Fix Dependency error

* Fix

* Fix index.d.ts

* Fix Documentation

* Fix Lint

* Fix index.d.ts

* Fix The Matrix parameter support and index.d.ts

* Fix Format

* Add Unit Tests for Matrix cases

* Fix Docs and Lint

* Add tests for coverage

* Fix Format

* Update AUTHORS

---------

Co-authored-by: Jos de Jong <wjosdejong@gmail.com>

AUTHORS

@@ -227,6 +227,7 @@ Maxim Mazurok <maxim@mazurok.com>
 kunalagrwl <kunalaggarwal82@gmail.com>
 Michael Greminger <michael.greminger@gmail.com>
 Kiku <kiku-cn@foxmail.com>
+Aly Khaled <alykhaled2001@live.com>
 MaybePixem <47889538+MaybePixem@users.noreply.github.com>
 
 # Generated by tools/update-authors.js

src/expression/embeddedDocs/embeddedDocs.js

@@ -181,6 +181,8 @@ import { setPowersetDocs } from './function/set/setPowerset.js'
 import { setSizeDocs } from './function/set/setSize.js'
 import { setSymDifferenceDocs } from './function/set/setSymDifference.js'
 import { setUnionDocs } from './function/set/setUnion.js'
+import { zpk2tfDocs } from './function/signal/zpk2tf.js'
+import { freqzDocs } from './function/signal/freqz.js'
 import { erfDocs } from './function/special/erf.js'
 import { madDocs } from './function/statistics/mad.js'
 import { maxDocs } from './function/statistics/max.js'
@@ -519,6 +521,10 @@ export const embeddedDocs = {
   setSymDifference: setSymDifferenceDocs,
   setUnion: setUnionDocs,
 
+  // functions - signal
+  zpk2tf: zpk2tfDocs,
+  freqz: freqzDocs,
+
   // functions - special
   erf: erfDocs,
 

src/expression/embeddedDocs/function/signal/freqz.js

@@ -0,0 +1,15 @@
+export const freqzDocs = {
+  name: 'freqz',
+  category: 'Signal',
+  syntax: [
+    'freqz(b, a)',
+    'freqz(b, a, w)'
+  ],
+  description: 'Calculates the frequency response of a filter given its numerator and denominator coefficients.',
+  examples: [
+    'freqz([1, 2], [1, 2, 3])',
+    'freqz([1, 2], [1, 2, 3], [0, 1])',
+    'freqz([1, 2], [1, 2, 3], 512)'
+  ],
+  seealso: []
+}

src/expression/embeddedDocs/function/signal/zpk2tf.js

@@ -0,0 +1,14 @@
+export const zpk2tfDocs = {
+  name: 'zpk2tf',
+  category: 'Signal',
+  syntax: [
+    'zpk2tf(z, p, k)'
+  ],
+  description: 'Compute the transfer function of a zero-pole-gain model.',
+  examples: [
+    'zpk2tf([1, 2], [-1, -2], 1)',
+    'zpk2tf([1, 2], [-1, -2])',
+    'zpk2tf([1 - 3i, 2 + 2i], [-1, -2])'
+  ],
+  seealso: []
+}

src/factoriesAny.js

@@ -259,6 +259,8 @@ export { createResolve } from './function/algebra/resolve.js'
 export { createSymbolicEqual } from './function/algebra/symbolicEqual.js'
 export { createDerivative } from './function/algebra/derivative.js'
 export { createRationalize } from './function/algebra/rationalize.js'
+export { createZpk2tf } from './function/signal/zpk2tf.js'
+export { createFreqz } from './function/signal/freqz.js'
 export { createReviver } from './json/reviver.js'
 export { createReplacer } from './json/replacer.js'
 export {

src/function/signal/freqz.js

@@ -0,0 +1,109 @@
+import { factory } from '../../utils/factory.js'
+
+const name = 'freqz'
+
+const dependencies = [
+  'typed',
+  'add',
+  'multiply',
+  'Complex',
+  'divide',
+  'matrix'
+]
+
+export const createFreqz = /* #__PURE__ */ factory(name, dependencies, ({ typed, add, multiply, Complex, divide, matrix }) => {
+  /**
+     * Calculates the frequency response of a filter given its numerator and denominator coefficients.
+     *
+     * Syntax:
+     *    math.freqz(b, a)
+     *    math.freqz(b, a, w)
+     *
+     * Examples:
+     *   math.freqz([1, 2], [1, 2, 3], 4) // returns { h: [0.5 + 0i, 0.4768589245763655 + 0.2861153547458193i, 0.25000000000000006 + 0.75i, -0.770976571635189 + 0.4625859429811135i], w: [0, 0.7853981633974483, 1.5707963267948966, 2.356194490192345 ] }
+     *   math.freqz([1, 2], [1, 2, 3], [0, 1]) // returns { h: [0.5 + 0i, 0.45436781 + 0.38598051i], w: [0, 1] }
+     *
+     * See also:
+     *  zpk2tf
+     *
+     * @param {Array.<number>} b The numerator coefficients of the filter.
+     * @param {Array.<number>} a The denominator coefficients of the filter.
+     * @param {Array.<number>} [w] A vector of frequencies (in radians/sample) at which the frequency response is to be computed or the number of points to compute (if a number is not provided, the default is 512 points)
+     * @returns {Object} An object with two properties: h, a vector containing the complex frequency response, and w, a vector containing the normalized frequencies (in radians/sample) at which the response was computed.
+     *
+     *
+     */
+  return typed(name, {
+    'Array, Array': function (b, a) {
+      const w = createBins(512)
+      return _freqz(b, a, w)
+    },
+    'Array, Array, Array': function (b, a, w) {
+      return _freqz(b, a, w)
+    },
+    'Array, Array, number': function (b, a, w) {
+      if (w < 0) {
+        throw new Error('w must be a positive number')
+      }
+      const w2 = createBins(w)
+      return _freqz(b, a, w2)
+    },
+    'Matrix, Matrix': function (b, a) {
+      // console.log('here')
+      const _w = createBins(512)
+      const { w, h } = _freqz(b.valueOf(), a.valueOf(), _w)
+      return {
+        w: matrix(w),
+        h: matrix(h)
+      }
+    },
+    'Matrix, Matrix, Matrix': function (b, a, w) {
+      const { h } = _freqz(b.valueOf(), a.valueOf(), w.valueOf())
+      return {
+        h: matrix(h),
+        w: matrix(w)
+      }
+    },
+    'Matrix, Matrix, number': function (b, a, w) {
+      if (w < 0) {
+        throw new Error('w must be a positive number')
+      }
+      const _w = createBins(w)
+      const { h } = _freqz(b.valueOf(), a.valueOf(), _w)
+      return {
+        h: matrix(h),
+        w: matrix(_w)
+      }
+    }
+  })
+
+  function _freqz (b, a, w) {
+    const num = []
+    const den = []
+    for (let i = 0; i < w.length; i++) {
+      let sumNum = Complex(0, 0)
+      let sumDen = Complex(0, 0)
+      for (let j = 0; j < b.length; j++) {
+        sumNum = add(sumNum, multiply(b[j], Complex(Math.cos(-j * w[i]), Math.sin(-j * w[i]))))
+      }
+      for (let j = 0; j < a.length; j++) {
+        sumDen = add(sumDen, multiply(a[j], Complex(Math.cos(-j * w[i]), Math.sin(-j * w[i]))))
+      }
+      num.push(sumNum)
+      den.push(sumDen)
+    }
+    const h = []
+    for (let i = 0; i < num.length; i++) {
+      h.push(divide(num[i], den[i]))
+    }
+    return { h, w }
+  }
+
+  function createBins (n) {
+    const bins = []
+    for (let i = 0; i < n; i++) {
+      bins.push(i / n * Math.PI)
+    }
+    return bins
+  }
+})

src/function/signal/zpk2tf.js

@@ -0,0 +1,86 @@
+import { factory } from '../../utils/factory.js'
+
+const name = 'zpk2tf'
+
+const dependencies = [
+  'typed',
+  'add',
+  'multiply',
+  'Complex',
+  'number'
+]
+
+export const createZpk2tf = /* #__PURE__ */ factory(name, dependencies, ({ typed, add, multiply, Complex, number }) => {
+  /**
+     * Compute the transfer function of a zero-pole-gain model.
+     *
+     * Syntax:
+     *      math.zpk2tf(z, p, k)
+     *
+     * Examples:
+     *    math.zpk2tf([1, 2], [-1, -2], 1)    // returns [[1, -3, 2], [1, 3, 2]]
+     *
+     * See also:
+     *   freqz
+     *
+     * @param {Array} z Array of zeros values
+     * @param {Array} p Array of poles values
+     * @param {number} k Gain value
+     * @return {Array} Two dimensional array containing the numerator (first row) and denominator (second row) polynomials
+     *
+     */
+  return typed(name, {
+    'Array,Array,number': function (z, p, k) {
+      return _zpk2tf(z, p, k)
+    },
+    'Array,Array': function (z, p) {
+      return _zpk2tf(z, p, 1)
+    },
+    'Matrix,Matrix,number': function (z, p, k) {
+      return _zpk2tf(z.valueOf(), p.valueOf(), k)
+    },
+    'Matrix,Matrix': function (z, p) {
+      return _zpk2tf(z.valueOf(), p.valueOf(), 1)
+    }
+  })
+
+  function _zpk2tf (z, p, k) {
+    // if z is bignumber, convert it to number
+    if (z.some((el) => el.type === 'BigNumber')) {
+      z = z.map((el) => number(el))
+    }
+    // if p is bignumber, convert it to number
+    if (p.some((el) => el.type === 'BigNumber')) {
+      p = p.map((el) => number(el))
+    }
+    let num = [Complex(1, 0)]
+    let den = [Complex(1, 0)]
+    for (let i = 0; i < z.length; i++) {
+      let zero = z[i]
+      if (typeof zero === 'number') zero = Complex(zero, 0)
+      num = _multiply(num, [Complex(1, 0), Complex(-zero.re, -zero.im)])
+    }
+    for (let i = 0; i < p.length; i++) {
+      let pole = p[i]
+      if (typeof pole === 'number') pole = Complex(pole, 0)
+      den = _multiply(den, [Complex(1, 0), Complex(-pole.re, -pole.im)])
+    }
+    for (let i = 0; i < num.length; i++) {
+      num[i] = multiply(num[i], k)
+    }
+    return [num, den]
+  }
+
+  function _multiply (a, b) {
+    const c = []
+    for (let i = 0; i < a.length + b.length - 1; i++) {
+      c[i] = Complex(0, 0)
+      for (let j = 0; j < a.length; j++) {
+        if (i - j >= 0 && i - j < b.length) {
+          c[i] = add(c[i], multiply(a[j], b[i - j]))
+        }
+      }
+    }
+    return c
+  }
+})

test/unit-tests/function/signal/freqz.test.js

@@ -0,0 +1,86 @@
+import approx from '../../../../tools/approx.js'
+import assert from 'assert'
+import math from '../../../../src/defaultInstance.js'
+
+const freqz = math.freqz
+
+describe('freqz', function () {
+  it('should return the frequency response of a zero-pole-gain model given number as w parameter', function () {
+    approx.deepEqual(
+      freqz([math.complex(1, 0), math.complex(-1, -5)], [math.complex(1, 0), math.complex(5, 0), math.complex(6, 0)], 5), {
+        h: [math.complex(0, -0.4166666666666667),
+          math.complex(0.08934733914100848, -0.3891571989056336),
+          math.complex(0.19350242447606233, -0.43679084380970035),
+          math.complex(0.5068505396172553, -0.5776505671214675),
+          math.complex(1.5607298559987381, -0.26338443553329166)
+        ],
+        w: [0, 0.62831853, 1.25663706, 1.88495559, 2.51327412]
+      })
+    approx.deepEqual(freqz([math.complex(1, 0), math.complex(-3, 0), math.complex(2, 0)], [math.complex(1, 0), math.complex(3, 0), math.complex(2, 0)], 5), {
+      h: [math.complex(0, 0),
+        math.complex(-0.09275445814522237, -0.11835248219923325),
+        math.complex(-0.4432171093183538, -0.3495195356882487),
+        math.complex(-1.3911165095967133, -1.0970298445163509),
+        math.complex(-4.102237436136192, -5.234357386651877)
+      ],
+      w: [0, 0.62831853, 1.25663706, 1.88495559, 2.51327412]
+    })
+  })
+
+  it('should return the frequency response of a zero-pole-gain model when not given w parameter', function () {
+    const { h, w } = freqz([math.complex(1, 0), math.complex(-1, -5)], [math.complex(1, 0), math.complex(5, 0), math.complex(6, 0)])
+    approx.deepEqual(h.length, 512)
+    approx.deepEqual(w.length, 512)
+  })
+
+  it('should return the frequency response of a zero-pole-gain model given b and a as matrix and not given w parameter', function () {
+    const b = math.matrix([math.complex(1, 0), math.complex(-1, -5)])
+    const a = math.matrix([math.complex(1, 0), math.complex(5, 0), math.complex(6, 0)])
+    const { h, w } = freqz(b, a)
+    approx.deepEqual(h._size, [512])
+    approx.deepEqual(w._size, [512])
+  })
+
+  it('should return the frequency response of a zero-pole-gain model given array as w parameter', function () {
+    approx.deepEqual(
+      freqz([math.complex(1, 0), math.complex(-1, -5)], [math.complex(1, 0), math.complex(5, 0), math.complex(6, 0)], [0, 1, 2]), {
+        h: [math.complex(0, -0.4166666666666667),
+          math.complex(0.1419346, -0.4055241),
+          math.complex(0.62506469, -0.59840473)
+        ],
+        w: [0, 1, 2]
+      })
+  })
+
+  it('should return the frequency response of a zero-pole-gain model given matrix as b,a and w parameter', function () {
+    approx.deepEqual(
+      freqz(math.matrix([math.complex(1, 0), math.complex(-1, -5)]), math.matrix([math.complex(1, 0), math.complex(5, 0), math.complex(6, 0)]), math.matrix([0, 1, 2])), {
+        h: math.matrix([math.complex(0, -0.4166666666666667),
+          math.complex(0.1419346, -0.4055241),
+          math.complex(0.62506469, -0.59840473)
+        ]),
+        w: math.matrix([0, 1, 2])
+      })
+  })
+
+  it('should return the frequency response of a zero-pole-gain model given matrix as b,a and number as w parameter', function () {
+    approx.deepEqual(
+      freqz(math.matrix([math.complex(1, 0), math.complex(-1, -5)]), math.matrix([math.complex(1, 0), math.complex(5, 0), math.complex(6, 0)]), 5), {
+        h: math.matrix([math.complex(0, -0.4166666666666667),
+          math.complex(0.08934733914100848, -0.3891571989056336),
+          math.complex(0.19350242447606233, -0.43679084380970035),
+          math.complex(0.5068505396172553, -0.5776505671214675),
+          math.complex(1.5607298559987381, -0.26338443553329166)
+        ]),
+        w: math.matrix([0, 0.62831853, 1.25663706, 1.88495559, 2.51327412])
+      })
+  })
+
+  it('should error with negative number of points', function () {
+    assert.throws(function () { freqz([1, 2], [1, 2, 3], -1) }, /w must be a positive number/)
+  })
+
+  it('should error with negative number of points when given matrix as b,a and w parameter', function () {
+    assert.throws(function () { freqz(math.matrix([1, 2]), math.matrix([1, 2, 3]), -1) }, /w must be a positive number/)
+  })
+})