-
Notifications
You must be signed in to change notification settings - Fork 54
/
designHBF7.m
executable file
·263 lines (246 loc) · 7.58 KB
/
designHBF7.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
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
function [f1_saved,f2_saved,info]=designHBF(fp,delta,debug)
%function [f1,f2,info]=designHBF(fp=0.2,delta=1e-5,debug=0)
%Design a half-band filter which can be realized without general multipliers.
%The filter is a composition of a prototype and sub- filter.
%Input
% fp The normalized cutoff frequency of the filter. Due to the
% symmetry imposed by a HBF, the stopband begins at 0.5-fp.
% delta The absolute value of the deviation of the frequency response from
% the ideal values of 1 in the passband and 0 in the stopband.
%
%Output
% f1,f2 The coefficients of the prototype and sub-filters
% and their canonical-signed digit (csd) representation.
% info A vector containing the following data (only set when debug=1):
% complexity The number of additions per output sample.
% n1,n2 The length of the f1 and f2 vectors.
% sbr The achieved stob-band attenuation (dB).
% phi The scaling factor for the F2 filter.
% To Do: Clean up the code a bit more, esp. wrt the use of the struct. arrays.
% Use the phi variable to cut down on the number of adders in F2.
% Apply a simulated annealing/genetic optimization alg instead
% of the ad hoc one I have now.
%Handle the input arguments
parameters = ['fp ';'delta';'debug'];
defaults = [ 0.2 1e-5 0];
for i=1:length(defaults)
if i>nargin
eval([parameters(i,:) '=defaults(i);'])
elseif eval(['any(isnan(' parameters(i,:) ')) | isempty(' parameters(i,:) ')'])
eval([parameters(i,:) '=defaults(i);'])
end
end
%Try several different values for the fp1 parameter.
%The best values are usually around .04
%Surrender if 3 successive attempts yield progressively greater complexity.
lowest_complexity = Inf; prev_complexity = Inf;
for fp1 = [.03 .035 .025 .040 .020 .045 .015 .05]
failed = 0;
[f1 zetap phi] = designF1( delta, fp1 );
if zetap == 1 % designF1 failed
failed = 1;
if debug
fprintf(2,'designF1 failed at fp1=%f\n',fp1);
end
end
if ~failed
f2 = designF2( fp, zetap, phi );
n1 = length(f1); n2 = length(f2);
if n2 == 0 % designF2 failed
failed = 1;
if debug
fprintf(2,'designF2 failed when zetap=%f, phi=%f\n',zetap,phi);
end
end
end
if ~failed
% complexity(+ performance) = the number of two-input adders (+ sbr)
complexity = size([f1.csd],2) + (2*n1-1)*(n2+size([f2.csd],2)-1);
if debug
msg = sprintf('%d adders: n1=%d, n2=%d, (fp1=%.2f, zetap=%.3f, phi=%4.2f)', ...
complexity, n1, n2, fp1, zetap, phi );
else
msg = '';
end
[fresp pbr sbr] = frespHBF([], f1, f2, phi, fp, msg);
if pbr <= delta & sbr <= delta
complexity = complexity + sbr;
if complexity < prev_complexity
worse = 0;
if complexity < lowest_complexity
lowest_complexity = complexity;
f1_saved = f1; f2_saved = f2;
phi_saved = phi;
if debug
fprintf( 1, '%s\n', msg )
end
end
else
worse = worse + 1;
if worse > 2
break;
end
end
prev_complexity = complexity;
end % if pbr <= delta
end
end % for fp1
if isinf(lowest_complexity)
fprintf(1,'%s: Unable to meet the design requirements.\n', mfilename);
elseif debug
complexity = floor(lowest_complexity);
msg = sprintf( 'Final Design: %d adders', complexity);
[junk pbr sbr] = frespHBF([], f1_saved, f2_saved, phi_saved, fp, msg);
n1 = length(f1_saved); n2 = length(f2_saved);
fprintf(1,'%s (%d,%d,%.0fdB)\n', msg,n1,n2,dbv(sbr));
info = [ complexity n1 n2 dbv(sbr) phi_saved ];
end
return
function [f1_saved,zetap,phi] = designF1(delta, fp1)
% [f1 zetap phi] = designF1(delta, fp1) Design the F1 sub-filter
% of a Saramaki halfband filter. This function is called by designHBF.m.
%
% f1 A structure array containing the F1 filter coefficents and
% Their CSD representation.
% phi The scaling factor for the F2 filter (imbedded in the f1 coeffs.)
passband = exp(4*pi*j*linspace(0,fp1));
ok = 0;
for n1 = 1:2:7 % Odd values only
if n1 == 1
h = [0.5 0.5];
else
h = firpm(2*n1-1,[0 4*fp1 1 1],[1 1 0 0]);
if ~(abs(sum(h)-1) < 1e-3 ) % remez bug! Use firls instead
h = firls(2*n1-1,[0 4*fp1 1-1e-6 1],[1 1 0 0]);
end
end
fresp = abs( polyval(h,passband) );
if max( abs(fresp-1) ) <= delta
ok = 1;
break
end
end
if ~ok
zetap = 1; % Use this as an indication that the function failed.
return
end
% Transform h(n) to a chebyshev polynomial f1(n)
% Sum(f1(i)*cos(w)^n)|i=1:n1 + Sum(h(n1+i))*cos(n*w))|i=1:n1, n = 2*i-1;
w = pi*rand(1,n1);
cos_w = cos(w);
A = zeros(n1,length(w));
B = zeros(1,n1);
for i = 1:n1
n = 2*i-1;
A(i,:) = cos_w .^ n;
B = B + h(n1+i)* cos(n*w);
end
f1 = B/A;
% Matlab Ver. 5 change:
phivecb = [];
% Optimize the quantized version of f1 to maximize the stopband width
% ( = acos(zetap) )
zetap = 1;
testPoints = [0 logspace(-2,0,128)] - 1;
for nsd = 3:8
f1a = f1'; f1b = f1'; % First try the unperturbed filter.
for phia = 1 ./ [1 f1]
phia = phia / 2^nextpow2(phia); % keep phi in (0.5,1]
% Try a bunch of coefficients in the current neighborhood,
% shrinking the neighborhood once 10 successive trial values show no
% improvement. If 2 successive shrinkages do no good, try a higher nsd.
count = 0;
nohelp = 0;
neighborhood = .05;
while neighborhood > 1e-5
phivec = phia .^ [1:2:2*n1-1]';
% Matlab Ver. 5 change:
if isempty(phivecb); phivecb = phivec; end
f1q = bquantize( f1a.*phivec, nsd );
F1 = evalF1( [f1q.val], testPoints, phia );
fi = find( abs(F1) > delta );
zeta = -testPoints( max( fi(1)-1, 1 ) );
%fprintf(2,'nsd=%d, nbhd= %f, count=%d, zeta = %f, phia=%f\n', ...
% nsd, neighborhood, count, zeta, phia );
if zeta < zetap
count = 0;
nohelp = 0;
zetap = zeta;
f1b = [f1q.val]';
f1_saved = f1q;
phi = phia;
phivecb = phivec;
else
count = count + 1;
end
if count > 10
count = 0;
neighborhood = neighborhood/2;
nohelp = nohelp +1;
if nohelp > 2
break;
end
end
f1a = f1b./phivecb + neighborhood*(rand(size(f1b))-0.5);
phia = phia + neighborhood*(rand(1,1)-0.5);
end
if zetap < 1 % Found a filter with adequate attn.
break;
end
end % for phia ...
if zetap < 1 % Found a filter with adequate attn.
break;
end
end
return
function f2 = designF2(fp,zetap,phi)
% f2 = designF2(fp,zetap,phi) Design the F2 sub-filter
% of a Saramaki halfband filter. This function is called by designHBF.m.
% subfilter design:
% 1 - delta2' < |F2/phi| < 1 for f in [0 fp];
% -1 < |F2/phi| < -1 + delta2' for f in [0.5-fp, 0.5];
% 1-delta2' = (1-delta2)/(1+delta2)
delta2 = (1-zetap)/(1+zetap);
%delta2p = 1 - (1-delta2)/(1+delta2);
% determine the minimum order required by the filter
passband = exp(j*linspace(0,4*pi*fp));
for nsub = 3:2:17
h2 = firpm(nsub,[0 4*fp 1 1], [1 1 0 0]);
mag = abs( polyval(h2,passband) );
if max(abs(mag-1)) < delta2;
break;
end
end
n2min = (nsub+1)/2;
% Search all n2,nsd pairs, in order of the product n2*(nsd+1)
% allowing fp to be a variable?
success = 0;
nsdmin = 3; nsdmax = 6;
for product = (nsdmin+1)*n2min:(nsdmax+1)*n2min
for nsd = nsdmin:nsdmax
n2 = product/(nsd+1);
if floor(n2) ~= n2 % Only take integer n2,nsd pairs
break
end
nsub = 2*n2-1;
% Could try a bunch of fp values
%fprintf(2,'designF2: Trying (n2,nsd2,fp)=(%2d,%2d,%6.4f)\n',n2,nsd,fp);
h2 = firpm(nsub,[0 4*fp 1 1], [1 1 0 0]);
h2 = h2/(phi*(1+delta2)); % Adjust the coefficients.
f2 = bquantize( h2(n2+1:nsub+1), nsd );
h2 = (1+delta2)*phi*[f2(n2:-1:1).val f2.val];
mag = abs( polyval(h2,passband) );
if max(abs(mag-1)) < delta2;
success =1;
break;
end
end
if success
break;
end
end
if ~success
f2 = [];
q2 = [];
end
return