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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
// 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/signal/
// Modifieded by: Abinash Singh , SOE CUSAT
// Last Modified on : Feb 2024
// Organization: FOSSEE, IIT Bombay
// Email: toolbox@scilab.in
/*
Calling Sequence
[b_out, a_out] = impinvar (b, a, fs, tol)
[b_out, a_out] = impinvar (b, a, fs)
[b_out, a_out] = impinvar (b, a)
Converts analog filter with coefficients b and a to digital, conserving impulse response.
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%) This function does the inverse of impinvar so that the following example should restore the original values of a and b.
invimpinvar implements the reverse of this function.
[b, a] = impinvar (b, a);
[b, a] = invimpinvar (b, a);
Reference: Thomas J. Cavicchi (1996) “Impulse invariance and multiple-order poles”. IEEE transactions on signal processing, Vol 44 (9): 2344–2347
Dependencies
residue
*/
function [b_out, a_out] = impinvar (b_in, a_in, fs , tol)
error("impinvar: Missing functionality .This function is not implemented yet. Will available in next release");
endfunction
/*
if (nargin <2)
error ("impinvar: Insufficient input arguments");
end
if nargin < 3
fs = 1;
end
if nargin < 4
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
[r_in, p_in, k_in] = residue(b_in, a_in); // partial fraction expansion
n = length(r_in); // Number of poles/residues
if (length(k_in)>0) // Greater than zero means we cannot do impulse invariance
error("Order numerator >= order denominator");
end
r_out = zeros(1,n); // Residues of H(z)
p_out = zeros(1,n); // Poles of H(z)
k_out = 0; // Constant term of H(z)
i=1;
while (i<=n)
m = 1;
first_pole = p_in(i); // Pole in the s-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, p, k] = z_res(r_in(i-m+1:i), first_pole, ts); // Find z-domain residues
k_out = k_out + k; // Add direct term to output
p_out(i-m+1:i) = p; // Copy z-domain pole(s) to output
r_out(i-m+1:i) = r; // Copy z-domain residue(s) to output
i=i+1; // Next s-domain residue/pole
end
[b_out, a_out] = inv_residue(r_out, p_out, k_out, tol);
a_out = to_real(a_out); // Get rid of spurious imaginary part
b_out = to_real(b_out);
// Shift results right to account for calculating in z instead of z^-1
b_out($)=[];
endfunction
*/
// Convert residue vector for single and multiple poles in s-domain (located at sm) to
// residue vector in z-domain. The variable k is the direct term of the result.
function [r_out, p_out, k_out] = z_res (r_in, sm, ts)
p_out = exp(ts * sm); // z-domain pole
n = length(r_in); // Multiplicity of the pole
r_out = zeros(1,n); // Residue vector
// First pole (no multiplicity)
k_out = r_in(1) * ts; // PFE of z/(z-p) = p/(z-p)+1; direct part
r_out(1) = r_in(1) * ts * p_out; // pole part of PFE
for i=(2:n) // Go through s-domain residues for multiple pole
r_out(1:i) = r_out(1:i) + r_in(i) * polyrev(h1_z_deriv(i-1, p_out, ts)); // Add z-domain residues
end
endfunction
function p_out = polyrev (p_in)
p_out = p_in($:-1:1);
endfunction
function p_out = to_real(p_in)
p_out = abs(p_in) .* sign(real(p_in));
endfunction
/*
//test passed
[b_out,a_out]=impinvar([1],[1 1],100);
assert_checkalmostequal(b_out,0.01,%eps,1e-4)
assert_checkalmostequal(a_out,[1 -0.99],%eps,1e-4)
//test passed
[b_out,a_out]=impinvar([1],[1 2 1],100)
assert_checkalmostequal(b_out,[0 9.9005e-5],%eps,1e-4)
assert_checkalmostequal(a_out,[1 -1.9801 0.9802],%eps,1e-4)
[b_out, a_out] = impinvar([1], [1 3 3 1], 100) // test passed
[b_out, a_out] = impinvar([1 1], [1 3 3 1], 100) // test passed
[b_out, a_out] = impinvar([1 1 1], [1 3 3 1], 100) // test passed
// FIXME : builtin filter doesn't accepts complex parameters
// These test cases will through errors
// [b_out, a_out] = impinvar([1], [1 0 1], 100)
// [b_out, a_out] = impinvar([1 1], [1 0 1], 100)
// [b_out, a_out] = impinvar([1], [1 0 2 0 1], 100)
// [b_out, a_out] = impinvar([1 1], [1 0 2 0 1], 100)
// [b_out, a_out] = impinvar([1 1 1], [1 0 2 0 1], 100)
// [b_out, a_out] = impinvar([1 1 1 1], [1 0 2 0 1], 100)
*/
|
