1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
// Copyright (C) 2018 - IIT Bombay - FOSSEE
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
// Original Source : https://octave.sourceforge.io/
// Modifieded by: Abinash Singh Under FOSSEE Internship
// Last Modified on : 3 Feb 2024
// Organization: FOSSEE, IIT Bombay
// Email: toolbox@scilab.in
// Impulse invariant conversion from s to z domain
/*
[b_out, a_out] = invimpinvar (b, a, fs, tol)
[b_out, a_out] = invimpinvar (b, a, fs)
[b_out, a_out] = invimpinvar (b, a)
Converts digital filter with coefficients b and a to analog, conserving impulse response.
This function does the inverse of impinvar so that the following example should restore the original values of a and b.
[b, a] = impinvar (b, a);
[b, a] = invimpinvar (b, a);
If fs is not specified, or is an empty vector, it defaults to 1Hz.
If tol is not specified, it defaults to 0.0001 (0.1%)
Dependencies
residue
inv_residue
*/
function [b_out, a_out] = invimpinvar (b_in, a_in, fs, tol)
error("invimpinvar: Missing functionality not implemented in this release .Will be Available soon ");
endfunction
// FIXME : fix filter function first . till then drop this fun
/*
if (nargin <2)
error("invimpinvar: Insufficient input arguments");
end
if nargin < 3 then fs = 1; end
if nargin < 4 then tol = 0.0001; end
// to be compatible with the matlab implementation where an empty vector can
// be used to get the default
if (isempty(fs))
ts = 1;
else
ts = 1/fs; // we should be using sampling frequencies to be compatible with Matlab
end
b_in = [b_in 0]; // so we can calculate in z instead of z^-1
[r_in, p_in, k_in] = residue(b_in, a_in); // partial fraction expansion
// clean r_in for zero values
n = length(r_in); // Number of poles/residues
if (length(k_in) > 1) // Greater than one means we cannot do impulse invariance
error("Order numerator > order denominator");
end
r_out = zeros(1,n); // Residues of H(s)
sm_out = zeros(1,n); // Poles of H(s)
i=1;
while (i<=n)
m=1;
first_pole = p_in(i); // Pole in the z-domain
while (i<n && abs(first_pole-p_in(i+1))<tol) // Multiple poles at p(i)
i=i+1; // Next residue
m=m+1; // Next multiplicity
end
[r, sm, k]= inv_z_res(r_in(i-m+1:i), first_pole, ts); // Find s-domain residues
k_in = k_in - k; // Just to check, should end up zero for physical system
sm_out(i-m+1:i) = sm; // Copy s-domain pole(s) to output
r_out(i-m+1:i) = r; // Copy s-domain residue(s) to output
i=i+1; // Next z-domain residue/pole
end
[b_out, a_out] = inv_residue(r_out, sm_out , 0, tol);
a_out = to_real(a_out); // Get rid of spurious imaginary part
b_out = to_real(b_out);
b_out = polyreduce(b_out);
endfunction
*/
// Inverse function of z_res (see impinvar source)
function [r_out, sm_out, k_out] = inv_z_res (r_in,p_in,ts)
n = length(r_in); // multiplicity of the pole
r_in = r_in.'; // From column vector to row vector
j=n;
while (j>1) // Go through residues starting from highest order down
r_out(j) = r_in(j) / ((ts * p_in)^j); // Back to binomial coefficient for highest order (always 1)
r_in(1:j) = r_in(1:j) - r_out(j) * polyrev(h1_z_deriv(j-1,p_in,ts)); // Subtract highest order result, leaving r_in(j) zero
j=j-1;
end
// Single pole (no multiplicity)
r_out(1) = r_in(1) / ((ts * p_in));
k_out = r_in(1) / p_in;
sm_out = log(p_in) / ts;
endfunction
/*
// tests passed
[b_out,a_out]=invimpinvar([1],[1 -0.5],0.01)
[b_out,a_out]=invimpinvar([1],[1 -1 0.25],0.01)
[b_out,a_out]=invimpinvar([1 1],[1 -1 0.25],0.01)
[b_out,a_out]=invimpinvar([1],[1 -1.5 0.75 -0.125],0.01)
[b_out,a_out]=invimpinvar([1 1],[1 -1.5 0.75 -0.125],0.01)
// FIXME : built in filter doesn't support complex parameters
// Because of this thsese test cases are failing
//[b_out,a_out]=invimpinvar([1],[1 0 0.25],0.01)
// [b_out,a_out]=invimpinvar([1 1],[1 0 0.25],0.01)
// [b_out,a_out]=invimpinvar([1],[1 0 0.5 0 0.0625],0.01)
// [b_out,a_out]=invimpinvar([1 1],[1 0 0.5 0 0.0625],0.01)
// [b_out,a_out]=invimpinvar([1 1 1],[1 0 0.5 0 0.0625],0.01
// [b_out,a_out]=invimpinvar([1 1 1 1],[1 0 0.5 0 0.0625],0.01)
*/
|
