dfildemo.m

数字通信第四版原书的例程

%DFILDEMO Nonlinear filter design problem using discrete optimization.
clf
echo on
% Example script-file for the design of finite precision filters
nbits = 8; 		% How many bits have we to realize filter
maxbin = 2^nbits-1; 	% Maximum number
n=4;			% Number of coefficients (Order of filter plus 1)
Wn = 1.5;		% Cut-off frequency for filter
Rp = 0.2;		% Decibels of ripple in the passband
w = 128; 		% Number of frequency points to take 

% Continuous filter design 
% (could use any of cheby, ellip, yulewalk or remez here.)
[b1,a1]=cheby1(n-1,Wn, Rp); 


pause % Strike any key to continue

[h,w]=freqz(b1,a1,w);	% Frequency response
h=abs(h);		% Magnitude response
plot(w, h), title('Frequency response using non-integer variables')
x = [b1,a1];		% The design variables

% Set bounds on the maximum and minimum values. 
if (any(x < 0))
% If there are negative coefficients - must use sign bit
	maxbin = floor(maxbin/2);
	vlb = -maxbin * ones(1, 2*n);
	vub = maxbin * ones(1, 2*n); 
else
% otherwise, all positive
	vlb = zeros(1,2*n); 
	vub = maxbin * ones(1, 2*n); 
end

% Set the biggest value equal to maxbin. 
[m, mix] = max(abs(x)); 
factor =  maxbin/m; 
x =  factor * x; 	% Rescale other filter coefficients
xorig = x;

xmask = 1:2*n;
% Remove the biggest value and the element that controls D.C. Gain
% from the list of values that can be changed. 
xmask([mix]) = [];
nx = 2*n; 

pause % Strike any key to continue

% Set termination criteria to reasonably high values 
% to promote fast convergence.

options = 1; 
options(2) = 0.1;  
options(3) = 1e-4;
options(4) = 1e-6; 

% Need to minimize absolute maximum values so set options(15) to the 
% number of frequency points. 

if length(w) == 1
	options(15) = w; 
else
	options(15) = length(w); 
end

pause % Strike any key to continue

% Discretize and eliminate first value
[x, xmask] = elimone(x, xmask, h, w, n, maxbin)

pause % Strike any key to continue

niters = length(xmask); 

pause % Strike any key to continue

disp(sprintf('Performing %g stages of optimization.\n\n', niters));

% This may take a long time ! (Ctrl-C to abort)

for m = 1:niters
	disp(sprintf('Stage: %g \n', m));
% Least squares solution:
%	x(xmask) = constr('filtfun2', x(xmask), options, vlb(xmask), vub(xmask), [], x, xmask, n, h, maxbin);


	x(xmask) = minimax('filtfun', x(xmask), options, vlb(xmask), vub(xmask), [], x, xmask, n, h, maxbin);
	[x, xmask] = elimone(x, xmask, h, w, n, maxbin);
end

% Check nearest integer values for better filter
pause % Strike any key to continue

xold = x;
xmask = 1:2*n;
xmask([n+1, mix]) = [];
x = x + 0.5; 
for i = xmask
	[x, xmask] = elimone(x, xmask, h, w, n, maxbin);
end
xmask = 1:2*n;
xmask([n+1, mix]) = [];
x= x - 0.5;
for i = xmask
	[x, xmask] = elimone(x, xmask, h, w, n, maxbin);
end
if any(abs(x) > maxbin), x = xold; end

% Frequency response of filter
pause % Strike any key to continue
subplot(211)
bo = x(1:n); 
ao = x(n+1:2*n); 
h2 = abs(freqz(bo,ao,128));
plot(w,h,w,h2,'o')
title('Optimized filter versus original')

% Compare it to a filter where the coefficients are just rounded 
% up or down.
xround =round(xorig)
b = xround(1:n); 
a = xround(n+1:2*n); 
h3 = abs(freqz(b,a,128));
subplot(212)
plot(w,h,w,h3,'+')
title('Rounded filter versus original')
set(gcf,'NextPlot','replace')