/
synthesizeNTF0.m
executable file
·141 lines (130 loc) · 3.05 KB
/
synthesizeNTF0.m
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
function ntf = synthesizeNTF0(order,OSR,opt,H_inf,f0)
% Determine the zeros.
if f0~=0 % Bandpass design-- halve the order temporarily.
order = order/2;
dw = pi/(2*OSR);
else
dw = pi/OSR;
end
if length(opt)==1
if opt==0
z = zeros(order,1);
else
z = dw*ds_optzeros(order,opt);
if isempty(z)
return;
end
end
if f0~=0 % Bandpass design-- shift and replicate the zeros.
order = order*2;
z = z + 2*pi*f0;
ztmp = [ z'; -z' ];
z = ztmp(:);
end
z = exp(j*z);
else
z = opt(:);
end
ntf = zpk(z,zeros(1,order),1,1);
itn_limit = 100;
% Iteratively determine the poles by finding the value of the x-parameter
% which results in the desired H_inf.
if f0 == 0 % Lowpass design
HinfLimit = 2^order; % !!! The limit is actually lower for opt=1 and low OSR
if H_inf >= HinfLimit
fprintf(2,'%s warning: Unable to achieve specified Hinf.\n', mfilename);
fprintf(2,'Setting all NTF poles to zero.\n');
ntf.p = zeros(order,1);
else
x=0.3^(order-1); % starting guess
converged = 0;
for itn=1:itn_limit
me2 = -0.5*(x^(2./order));
w = (2*[1:order]'-1)*pi/order;
mb2 = 1+me2*exp(j*w);
p = mb2 - sqrt(mb2.^2-1);
out = find(abs(p)>1);
p(out) = 1./p(out); % reflect poles to be inside the unit circle.
ntf.z = z; ntf.p = cplxpair(p);
f = real(evalTF(ntf,-1))-H_inf;
% [ x f ]
if itn==1
delta_x = -f/100;
else
delta_x = -f*delta_x/(f-fprev);
end
xplus = x+delta_x;
if xplus>0
x = xplus;
else
x = x*0.1;
end
fprev = f;
if abs(f)<1e-10 | abs(delta_x)<1e-10
converged = 1;
break;
end
if x>1e6
fprintf(2,'%s warning: Unable to achieve specified Hinf.\n', mfilename);
fprintf(2,'Setting all NTF poles to zero.\n');
ntf.z = z; ntf.p = zeros(order,1);
break;
end
if itn == itn_limit
fprintf(2,'%s warning: Danger! Iteration limit exceeded.\n',...
mfilename);
end
end
end
else % Bandpass design.
x = 0.3^(order/2-1); % starting guess (not very good for f0~0)
if f0>0.25
z_inf=1;
else
z_inf=-1;
end
c2pif0 = cos(2*pi*f0);
for itn=1:itn_limit
e2 = 0.5*x^(2./order);
w = (2*[1:order]'-1)*pi/order;
mb2 = c2pif0 + e2*exp(j*w);
p = mb2 - sqrt(mb2.^2-1);
% reflect poles to be inside the unit circle.
out = find(abs(p)>1);
p(out) = 1./p(out);
ntf.z = z; ntf.p = cplxpair(p);
f = real(evalTF(ntf,z_inf))-H_inf;
% [x f]
if itn==1
delta_x = -f/100;
else
delta_x = -f*delta_x/(f-fprev);
end
xplus = x+delta_x;
if xplus > 0
x = xplus;
else
x = x*0.1;
end
fprev = f;
if abs(f)<1e-10 | abs(delta_x)<1e-10
break;
end
if x>1e6
fprintf(2,'%s warning: Unable to achieve specified Hinf.\n', mfilename);
fprintf(2,'Setting all NTF poles to zero.\n');
ntf.p = zeros(order,1);
break;
end
if itn == itn_limit
fprintf(2,'%s warning: Danger! Iteration limit exceeded.\n',...
mfilename);
end
end
end
z = cplxpair(ntf.z{:});
p = cplxpair(ntf.p{:});
rev = order:-1:1;
% Assemble the ntf struct
ntf.z = z(rev)';
ntf.p = p(rev)';