home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 4.ddi / OPTIM.DI$ / DFILDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.3 KB  |  129 lines
  1. %DFILDEMO Nonlinear filter design problem using discrete optimization.
  2. clf
  3. echo on
  4. % Example script-file for the design of finite precision filters
  5. nbits = 8;         % How many bits have we to realize filter
  6. maxbin = 2^nbits-1;     % Maximum number
  7. n=4;            % Number of coefficients (Order of filter plus 1)
  8. Wn = 1.5;        % Cut-off frequency for filter
  9. Rp = 0.2;        % Decibels of ripple in the passband
  10. w = 128;         % Number of frequency points to take 
  11.  
  12. % Continuous filter design 
  13. % (could use any of cheby, ellip, yulewalk or remez here.)
  14. [b1,a1]=cheby1(n-1,Wn, Rp); 
  15.  
  16.  
  17. pause % Strike any key to continue
  18.  
  19. [h,w]=freqz(b1,a1,w);    % Frequency response
  20. h=abs(h);        % Magnitude response
  21. plot(w, h), title('Frequency response using non-integer variables')
  22. x = [b1,a1];        % The design variables
  23.  
  24. % Set bounds on the maximum and minimum values. 
  25. if (any(x < 0))
  26. % If there are negative coefficients - must use sign bit
  27.     maxbin = floor(maxbin/2);
  28.     vlb = -maxbin * ones(1, 2*n);
  29.     vub = maxbin * ones(1, 2*n); 
  30. else
  31. % otherwise, all positive
  32.     vlb = zeros(1,2*n); 
  33.     vub = maxbin * ones(1, 2*n); 
  34. end
  35.  
  36. % Set the biggest value equal to maxbin. 
  37. [m, mix] = max(abs(x)); 
  38. factor =  maxbin/m; 
  39. x =  factor * x;     % Rescale other filter coefficients
  40. xorig = x;
  41.  
  42. xmask = 1:2*n;
  43. % Remove the biggest value and the element that controls D.C. Gain
  44. % from the list of values that can be changed. 
  45. xmask([mix]) = [];
  46. nx = 2*n; 
  47.  
  48. pause % Strike any key to continue
  49.  
  50. % Set termination criteria to reasonably high values 
  51. % to promote fast convergence.
  52.  
  53. options = 1; 
  54. options(2) = 0.1;  
  55. options(3) = 1e-4;
  56. options(4) = 1e-6; 
  57.  
  58. % Need to minimize absolute maximum values so set options(15) to the 
  59. % number of frequency points. 
  60.  
  61. if length(w) == 1
  62.     options(15) = w; 
  63. else
  64.     options(15) = length(w); 
  65. end
  66.  
  67. pause % Strike any key to continue
  68.  
  69. % Discretize and eliminate first value
  70. [x, xmask] = elimone(x, xmask, h, w, n, maxbin)
  71.  
  72. pause % Strike any key to continue
  73.  
  74. niters = length(xmask); 
  75.  
  76. pause % Strike any key to continue
  77.  
  78. disp(sprintf('Performing %g stages of optimization.\n\n', niters));
  79.  
  80. % This may take a long time ! (Ctrl-C to abort)
  81.  
  82. for m = 1:niters
  83.     disp(sprintf('Stage: %g \n', m));
  84. % Least squares solution:
  85. %    x(xmask) = constr('filtfun2', x(xmask), options, vlb(xmask), vub(xmask), [], x, xmask, n, h, maxbin);
  86.  
  87.  
  88.     x(xmask) = minimax('filtfun', x(xmask), options, vlb(xmask), vub(xmask), [], x, xmask, n, h, maxbin);
  89.     [x, xmask] = elimone(x, xmask, h, w, n, maxbin);
  90. end
  91.  
  92. % Check nearest integer values for better filter
  93. pause % Strike any key to continue
  94.  
  95. xold = x;
  96. xmask = 1:2*n;
  97. xmask([n+1, mix]) = [];
  98. x = x + 0.5; 
  99. for i = xmask
  100.     [x, xmask] = elimone(x, xmask, h, w, n, maxbin);
  101. end
  102. xmask = 1:2*n;
  103. xmask([n+1, mix]) = [];
  104. x= x - 0.5;
  105. for i = xmask
  106.     [x, xmask] = elimone(x, xmask, h, w, n, maxbin);
  107. end
  108. if any(abs(x) > maxbin), x = xold; end
  109.  
  110. % Frequency response of filter
  111. pause % Strike any key to continue
  112. subplot(211)
  113. bo = x(1:n); 
  114. ao = x(n+1:2*n); 
  115. h2 = abs(freqz(bo,ao,128));
  116. plot(w,h,w,h2,'o')
  117. title('Optimized filter versus original')
  118.  
  119. % Compare it to a filter where the coefficients are just rounded 
  120. % up or down.
  121. xround =round(xorig)
  122. b = xround(1:n); 
  123. a = xround(n+1:2*n); 
  124. h3 = abs(freqz(b,a,128));
  125. subplot(212)
  126. plot(w,h,w,h3,'+')
  127. title('Rounded filter versus original')
  128. set(gcf,'NextPlot','replace')
  129.