b2=firquest(n,ns,Fp,Fst,Ast,W,N)

% FIRQUEST Generate lowpass FIR filter with quantized coefficients
% Usage:
%   b = FIRQUEST(n,ns,Fp,Fst,Ast,W,N) 
% Arguments:
%   n   Lowpass Filter order 
%   ns  Shaping Filter order
%   Fp  passband frequency (normalized to Nyquist)
%   Fst stopband frequency (normalized to Nyquist)
%   Ast stopband suppression (dB)
%   W   coefficient word size in bits
%   N   fft size (defaults to 1024)
% Return values:
%   b   quantized filter coefficients
% Example:
%   b=firquest(37,6,0.3,0.6,16,120,1024);
%   fvtool(b,1);
% Based on:
%   J.J. Nielsen, 'Design of Linear-Phase Direct-Form FIR Digital Filters 
%   with Quantized Coefficients Using Error Spectrum Shaping Techniques', 
%   IEEE Transactions on Acoustics, Speech and Signal Processing, 
%   vol. ASSP-37, pp. 1020-1026, july 1989.

function b2=firquest(n,ns,Fp,Fst,Ast,W,N)

% default values
if ~exist('n'  ,'var');  n=37;       end
if ~exist('ns' ,'var');  ns=6;       end
if ~exist('Fp' ,'var');  Fp=0.3;     end
if ~exist('Fst','var');  Fst=0.6;    end
if ~exist('q'  ,'var');  W=16;       end
if ~exist('Ast','var');  Ast=W*6+40; end
if ~exist('N'  ,'var');  N=1024;     end

% design the lowpass filter
fd=fdesign.lowpass('n,fp,fst,ast',n,Fp,Fst,Ast);
d=design(fd,'equiripple');
b0=d.Numerator;

% design the shaping filter
dp=1E-1;
ds=1E-2;
bb=fircls1(ns,(Fp+Fst)/2,dp,ds,Fp,Fst,1E4);

%fd=fdesign.lowpass('n,fp,fst',3,Fp,Fst);
%d=design(fd,'equiripple');
%bb=d.Numerator;

% determine the stopband indices
f=linspace(0,2*(N-1)/N,N);
ist=find(f>Fst & f<1);

% rescaling
r=2^(W-1);
b0=r*b0/max(abs(b0))*(r-10)/r;
s=sum(b0);

% normal rounding
b1=round(b0);
y1=firquest_eval(b1,N,ist,s);
b2=b1;
y2=y1;
for i=1:10000
    bt=firquest_shape(bb,b0);
    y=firquest_eval(bt,N,ist,s);
    if y<y2 
        y2=y;
        b2=bt;
        display(sprintf('%5i %7.2f %5.2f',i,y2,y1-y2));
    end
end


if nargout==0
  h_bb=freqz(bb/bb(1),1,pi*f);
  h_b0=freqz(b0/s    ,1,pi*f);
  h_b1=freqz(b1/s    ,1,pi*f);
  h_b2=freqz(b2/s    ,1,pi*f);

  figure(1);
  plot(f,20*log10(abs([h_b1;h_b2;h_b0])));
  axis([0 1 -200 50]);
  grid on;
  title('Magnitude Response');
  xlabel('Normalized Frequency');
  ylabel('Magnitude (dB)');
  legend( 'round', 'shaped', 'unquantized');

  figure(2)
  plot(f,20*log10(abs([h_b1-h_b0;h_b2-h_b0;h_bb])));
  axis([0 1 -200 50]);
  grid on;
  title('Error Response');
  xlabel('Normalized Frequency');
  ylabel('Error (dB)');
  legend('round', 'shaped', 'shaping filter');
  clear b2
end

% evaluate filter performance by finding the max stopband level
function y=firquest_eval(bt,N,ist,s)
bf=fft([bt zeros(1,N-length(bt))]);
bst=bf(ist);
y=10*log10(max(bst.*conj(bst))/(s*s));

% generate a shaped filter, using the error spectral shaping technique
function z=firquest_shape(b,x)
L=length(x);
% even and odd lengths are treated seperately
if mod(L,2)==0
    x=x((L+2)/2:L);
else
    x=x((L+1)/2:L);
end

b=b/b(1);
b=b(2:length(b));
y=(rand(size(b))-0.5);
z=zeros(size(b));
for i=1:length(x)
    y0=sum(b.*y)+x(i);
    yr=round(y0);    
    z(i)=yr;
    y=[yr-y0 y(1:length(y)-1)];
end
if mod(L,2)==0 
    z=[fliplr(z) z];
else
    z=[fliplr(z) z(2:length(z))];
end