View Raw SPL
/*****************************************************************************
* *
* BILINEAR.SPL Copyright (C) 1998-2009 DSP Development Corporation *
* All Rights Reserved *
* *
* Author: Matthew Tom *
* *
* Synopsis: Bilinear transformation with optional frequency prewarping. *
* *
* Revisions: 3 May 1997 MAT Creation *
* 4 Jun 1998 MAT Restructuring *
* 8 Jun 1998 MAT Help Menu Entry *
* 23 Nov 2009 RRR Z transform form *
* *
*****************************************************************************/
#include
#if @HELP_BILINEAR
BILINEAR
Purpose: Converts a continuous system to a discrete system.
Syntax: (Zd,Pd,Kd) = Bilinear(Z,P,K,Fs,Fp)
Z - A row vector with the zeros of the transfer function.
P - A row vector with the poles of the transfer function.
K - A scalar gain constant.
Fs - A real specifying the sampling frequency in Hertz
Fp - Optional. Frequency in Hertz specifying at which point
the frequency responses before and after mapping match
exactly.
State-Space Syntax:
(Ad,Bd,Cd,Dd) = Bilinear(A,B,C,D,Fs,Fp)
A,B,C,D - Arrays containing the state-space representation of
the filter to be transformed.
Fs - A real specifying the sampling frequency in Hertz
Fp - Optional. Frequency in Hertz specifying at which point
the frequency responses before and after mapping match
exactly.
Transfer Function Syntax:
(NUMd,DENd) = Bilinear(NUM,DEN,Fs,Fp)
NUM - Column vector containing the numerator transfer
function coefficients in descending powers of s.
DEN - Column vector containing the denominator transfer
function coefficients in descending powers of s.
Fs - A real specifying the sampling frequency in Hertz
Fp - Optional. Frequency in Hertz specifying at which point
the frequency responses before and after mapping match
exactly.
Returns:
For the Zero-Pole-Gain format, two series and a constant.
Zd is the vector of zeros, Pd is the series of poles, and
Kd is the gain constant.
For the State-Space format, Ad, Bd, Cd and Dd are the
state-space representations of the z-transform discrete
equivalent filter.
For the Transfer Function format, B and A are series
containing numerator and denominator transfer function
coefficients.
Example:
b = {100};
a = {1, 2, 100};
(b1, a1) = bilinear(b, a, 200);
W1: (h, f) = sfreq(b, a, 1024);xy(f, mag(h))
W2: (h1, f1) = zfreq(b1, a1, 8192, 200);xy(f1, mag(h1))
W1 displays the magnitude of the frequency response for the
continuous system:
100
H(s) = -------------
2
s + 2s + 100
The system is converted to a digital system with a sample
rate of 200 Hz where:
b1 == {0.000622, 0.001243, 0.000622}
a1 == {1.000000, -1.987570, 0.990056}
representing the system:
-1 -2
0.000622 + 0.001243 z + 0.000622 z
H(z) = --------------------------------------
-1 -2
1 - 1.987570 z + 0.990056 z
W2 displays the magnitude of the frequency response for the
digital system.
Remarks:
The bilinear transform converts a continuous, constant
coefficient, time invariant analog system to a discrete,
constant coefficient shift invariant digital system. In
particular, the bilinear transform converts the analog
system Ha(s) to a digital system Hd(s) with a sample
period Ts = 1/Fs using the following mapping:
H(z) = H(s) |
| s = 2*Fs*(z-1)/(z+1)
For the Zero-Pole-Gain format, BILINEAR maps the zeros Z,
poles P and gain K of the continuous system to zeros Zd,
poles Pd and gain Kd the discrete system.
For the State-Space format, BILINEAR maps the continuous
state space matrices A, B, C, D where:
.
x = Ax + Bu
y = Cx + Du
to the discrete state space matrices Ad, Bd, Cd, Dd where:
x[n+1] = Ad x[n] + Bd u[n]
y[n] = Cd x[n] + Dd u[n]
For the Transfer Function format, BILINEAR maps the
continuous system H(s) to the discrete system H(z) where:
B[1] s^(N-1) + B[2] s^(N-2) + ... + B[N]
H(s) = ----------------------------------------
A[1] s^(M-1) + A[2] s^(M-2) + ... + A[M]
b[1] + b[2] z^(-1) + ... + b[N] z^(-N+1)
H(z) = ----------------------------------------
a[1] + a[2] z^(-1) + ... + a[M] z^(-M+1)
See Also:
FILTEQ
SFREQ
ZFREQ
ZPFCOEF
#endif
/* map H(s) to H(z) */
bilinear(z, p, k, fs, fp, fp1)
{
local mn, nn, md, nd, pd, zd, kd;
local a, b, c, d, n, t, r, t1, t2, ad, bd, cd, dd, num, den, NUMd, DENd;
(mn, nn) = size(z);
(md, nd) = size(p);
if (outargc == 3)
{
/* inputs in zero-pole-gain form */
if (mn > md)
{
error("bilinear - numerator order cannot exceed denominator");
}
if (argc == 5)
{
/* prewarp */
fp = 2 * pi * fp;
fs = fp / tan(fp / (2 * fs));
}
else
{
fs = 2 * fs;
}
/* remove nan and inf but preserve if empty */
if (not(isempty(z)))
{
z = delete(z, not(finite(z)));
}
pd = (1 + p / fs) / (1 - p / fs);
zd = (1 + z / fs) / (1 - z / fs);
kd = real(k * prod(fs - z) / prod(fs - p));
zd = {zd, -ones(length(pd) - length(zd), 1)};
return(zd, pd, kd);
}
else if (outargc == 4)
{
/* state-space case */
a = z;
b = p;
c = k;
d = fs;
fs = fp;
if (argc == 6)
{
/* prewarp */
fp = fp1;
fp = 2 * pi * fp;
fs = fp / (2 * tan(fp / (2 * fs)));
}
/* state-space version of the bilinear transformation */
n = max(size(a));
t = 1 / fs;
r = sqrt(t);
t1 = eye(n) + a * t / 2;
t2 = eye(n) - a * t / 2;
ad = t2 \^ t1;
bd = (t / r) * (t2 \^ b);
cd = r * (t2' \^ c')';
dd = (c' \^ t2')' *^ b * t / 2 + d;
return(ad, bd, cd, dd);
}
else if ((nd == 1) && (nn == 1))
{
if (mn > md)
{
error("bilinear - order of numerator cannot exceed order of denominator");
}
num = z;
den = p;
if (argc == 4)
{
fp = fs;
fs = k;
fp = 2 * pi * fp;
fs = fp / (2 * tan(fp / (2 * fs)));
}
else
{
fs = k;
}
/* convert to state-space form */
(a, b, c, d) = tf2ss(num, den);
/* state-space bilinear transformation */
n = max(size(a));
t = 1 / fs;
r = sqrt(t);
t1 = eye(n) + a * t / 2;
t2 = eye(n) - a * t / 2;
ad = t2 \^ t1;
bd = (t / r) * (t2 \^ b);
cd = r * (t2' \^ c')';
dd = (t2' \^ c')' *^ b * t / 2 + d;
/* convert to transfer function form */
DENd = poly(ad);
NUMd = poly(ad - bd *^ cd) + (dd[1] - 1) * DENd;
return(NUMd, DENd);
}
else
{
error("bilinear - first two arguments must be of the same dimension");
}
}