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
147
|
function L = filternorm(b,a,varargin)
// Calculates the L-2 norm or L-infinity norm of a digital filter
//
// Calling Sequence
// L = filternorm(b,a)
// L = filternorm(b,a,pnorm)
// L = filternorm(b,a,2,tol)
//
//
// Parameters
// b: The filter numerator coefficients.
// a: The filter denominator coefficients.
// pnorm: The L-norm to be calculated. The values accepted are 2 (L2 norm) or %inf (L-infinity norm). Default value is 2.
// tol: The tolerance of the L-2 norm to be calculated. If tol not specified, it defaults to 10^(-8). tol must be a positive scalar
//
//
// Examples
// // 1) L-2 norm of an IIR filter with tol = 10^(-10)
// b = [-3 2];
// a = [1 -0.5];
// L = filternorm(b, a, 2, 10d-10);
//
//
// See also
// norm
// zp2sos
//
// Authors
// Ayush Baid
exec('impz.sci', -1);
// ** Check on number of input, output arguments
[numOutArgs, numInArgs] = argn(0);
if numInArgs<2 | numInArgs>4 then
msg = "filternorm: Wrong number of input argument; 2-4 expected";
error(77,msg);
end
if numOutArgs~=1 then
msg = "filternorm: Wrong number of output argument; 1 expected";
error(78,msg);
end
// ** Check on b and a **
if isempty(a) then
a = 1;
end
if isempty(b) then
b = 1;
end
b = b(:);
a = a(:);
// check on datatype
if type(b)~=1 & type(b)~=8 then
msg = "filternorm: Wrong type for argument #1 (b): Real or complex matrix expected";
error(53,msg);
end
if type(a)~=1 & type(a)~=8 then
msg = "filternorm: Wrong type for argument #2 (a): Real or complex matrix expected";
error(53,msg);
end
// check on dimensions
if size(b,1)==1 then
b = b(:);
end
if size(a,1)==1 then
a = a(:);
end
if size(b,2)~=size(a,2) then
msg = "filternorm: Wrong size for arguments #1 (b) and #2(a): Same number of columns expected";
error(60,msg);
end
// ** Parsing the remaining arguments **
if length(varargin)==1 & ~isempty(varargin) then
pnorm = varargin(1);
tol = 1e-8;
elseif length(varargin)==2 then
pnorm = varargin(1);
tol = varargin(2);
if tol<=0 | length(tol)~=1 then
msg = "filternorm: Wrong value for argument #4 (tol): Must be a positive real scalar";
error(116,msg);
end
else
pnorm = 2;
tol = 1e-8;
end
if pnorm~=2 & length(varargin)==2 then
msg = "filternorm: Warning - Wrong value for argument #3 (pnorm): Must be 2 when tolerance is used";
end
// ** Calculations **
if isinf(pnorm) then
// We need to compute the frequency response and then get the one
// with the highest magnitude
h = frmag(b, a, 1024);
L = max(h);
else
if size(a,1) == 1 then
// the filter is FIR; impluse response is the filter coeffs
L = norm(b,pnorm)/a;
else
// the filter is IIR
// Checking for stability, as we wont be able to calc impulse response
// within a given tolerance.
pole_mag = abs(roots(a));
// stability check
max_dist = max(pole_mag);
if max_dist>=1 then
// poles lie on the unit circle or outside it. We do not have a
// decaying impulse response and hence truncation is not advisable
msg = "filternorm: Non convergent impulse response. All poles should lie inside the unit circle";
error(msg);
end
// ****
// Theory: (assuming stable filter)
// Each pole will contribute a decaying exponential. The pole with
// the highest magnitude will decay the slowest (i.e. will be the most
// dominant). Therefore, we will work with pole(s) having the largest
// magnitude to obtain a bound on the L2 norm of the tail.
// ****
// get the multiplicity of the largest pole
mult = sum(pole_mag>(max_dist-1e-5) & pole_mag<(max_dist+1e-5));
// Using integration of a^(-x) to get a bound
N = mult*log(tol)/log(max_dist);
// TODO: get filter coeffs using impzlength from octave
[h, temp1] = impz(b,a);
L = norm(h,2);
end
end
endfunction
|
