zpk2tf and freqz functions implementation

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/)
+ })
+})