-
Notifications
You must be signed in to change notification settings - Fork 50
/
synthesizeQNTF.m
executable file
·141 lines (139 loc) · 3.51 KB
/
synthesizeQNTF.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 = synthesizeQNTF(order,OSR,f0,NG,ING,n_im)
%ntf = synthesizeQNTF(order=3,OSR=64,f0=0,NG=-60,ING=-20,n_im=floor(order/3))
%Synthesize a noise transfer function for a quadrature delta-sigma modulator.
% order order of the modulator
% OSR oversampling ratio
% f0 center frequency (1->fs)
% NG in-band noise gain (dB)
% ING image-band noise gain (dB)
% n_im number of in-band image zeros
%
%ntf is a zpk object containing the zeros and poles of the NTF. See zpk.m
%
% ALPHA VERSION
% This function uses an experimental ad hoc method that is
% neither optimal nor robust.
% Handle the input arguments
parameters = {'order','OSR','f0','NG','ING','n_im'};
defaults = { 4 64 0 -60 -20 []};
for arg_i=1:length(defaults)
parameter = char(parameters(arg_i));
if arg_i>nargin | ( eval(['isnumeric(' parameter ') ']) & ...
eval(['any(isnan(' parameter ')) | isempty(' parameter ') ']) )
eval([parameter '=defaults{arg_i};'])
end
end
if isempty(n_im)
n_im = floor(order/3);
end
debug_it = 0;
if n_im==0
% Use synthesizeNTF to get an NTF with the specified NG; ignore ING
f1 = 0.5/OSR;
x = 1.5;
lowest_f = Inf;
dfdx = NaN;
for itn = 1:20
ntf = synthesizeNTF(order,OSR,1,x);
f = dbv(rmsGain(ntf,0,f1)) - NG;
if debug_it
fprintf(1,'x=%.2f f=%.2f\n', x, f);
end
if abs(f) < 0.01
break;
end
if isnan(dfdx)
dx = 0.1*sign(f);
dfdx = 0;
else
dfdx = (f-f_old)/dx;
dx_old = dx;
dx = -f/dfdx;
if abs(dx) > max(1,2*abs(dx_old))
dx = sign(dx)*max(1,2*abs(dx_old));
end
if x+dx <=1 % Hinf must be at least 1
dx = dx/2;
end
end
f_old = f;
x = x + dx;
end
if itn==20
fprintf(1,'Warning: Iteration limit reached. NTF may be poor\n');
end
% Rotate the NTF
z0 = exp(2i*pi*f0);
ntf = zpk(z0*ntf.z{1}, z0*ntf.p{1}, ntf.k, 1);
else
n_in = order-n_im;
f1 = f0-0.5/OSR;
f2 = f0+0.5/OSR;
z0 = exp(2i*pi*f0);
x = [20 20]; % "R" parameters for cheby2()
lowest_f = Inf;
dfdx = [NaN NaN];
for itn = 1:20
if debug_it
fprintf(1,'x=[%.2f %.2f], ', x(1), x(2) );
end
[b1 a1] = cheby2(n_in,x(1),1/OSR,'high');
[b2 a2] = cheby2(n_im,x(2),1/OSR,'high');
warning('off');
ntf0 = zpk( [roots(b1)*z0; roots(b2)*conj(z0)], ...
[roots(a1)*z0; roots(a2)*conj(z0)],1,1);
freq = linspace(-0.5,0.5,200);
m = evalTF(ntf0,exp(2i*pi*freq));
NG0 = dbv(rmsGain(ntf0,f1,f2));
ING0 = dbv(rmsGain(ntf0,-f1,-f2));
if debug_it
clf;
subplot(121);
plotPZ(ntf0);
subplot(122);
fprintf(1,'NG=%.1f, ING=%.1f\n', NG0, ING0);
plot(freq,dbv(m));
figureMagic([-0.5,0.5],0.05,2, [-100 30],10,2)
hold on;
plot([f1 f2],[1 1]*NG0,'k');
text(mean([f1 f2]), NG0, sprintf('NG=%.1fdB',NG0),'vert','bot');
plot([-f1 -f2],[1 1]*ING0,'k');
text(mean(-[f1 f2]), ING0, sprintf('ING=%.1fdB',ING0),'vert','bot');
drawnow;
end
f = [NG0 ING0] - [NG ING];
if abs(f) < 0.01
break;
end
if norm(f)<lowest_f
lowest_f = norm(f);
best = ntf0;
end
if abs(f(1)) > abs(f(2))
i=1; % adjust x(1)
else
i=2; % adjust x(2)
end
if isnan(dfdx(i))
dx = sign(f(i));
dfdx(i) = 0;
dfdx(3-i) = NaN;
else
dfdx(i) = (f(i)-f_old(i))/dx;
dfdx(3-i) = NaN;
dx = -f(i)/dfdx(i);
xnew = x(i) + dx;
if xnew < 0.5*x(i)
dx = -0.5*x(i);
elseif xnew > 2*x(i)
dx = x(i);
end
end
f_old = f;
x(i) = x(i) + dx;
end
if itn==20
fprintf(1,'Warning: Iteration limit reached. NTF may be poor\n');
end
ntf = best;
end