function  [h3]=comb112_sine(fi)

%      System decimates by 10
%      Function computes equalized frequency response of 8-fold
%      cascade of comb15 --> total comb filter order = 112
%      Recursive implementation requires shift by 31 bits at input
%      This yiels an efffective comb filter gain of G=1.1934
% m3 = number of fixed passband gain values for equalizer
% b0 = coefficient resolution of sharpening filter
% d1 = decimation factor
% fs = input sampling frequency in kHz
% fi = input frequency in kHz

b0=21;
b0=2^(b0-8);          % for -2^7 < coeff < +2^7
m3=13;                % optimum zero loss frequency spacing for m3=16
  z3=[.040 .122 .206 .285 .372 .449   .540 .617 .706 ...
      .784 .871 .940 .993];

n=2^16;               % fft length
fs=25000;	      % input sampling frequency in kHz
fp=1000;	      % filter passband edge frequency in kHz
d1=10;		      % decimation factor of comb filter
k1=15;                % order of comb filter
N1=8*(k1-1)           % order of cascaded comb filter
f1=(1:n)*(fs/(n*d1));
n2=floor(n/d1);
f2=(1:n2)*(fs/(n2*d1));

% Creat sinusoidal input
no1(1:n)=-.5+rand(n,1);   % noise with constant pdf between [-.5 .5]
wn=(1:n)*2*pi*fi/fs;
x=sin(wn)+no1*3e-5;       % sinusoidal input with additive noise
x=.7*x*(15/16)^8*2;       % gain corrected input due to scaling facor

% Compute cascade of comb filter impulse response
% 8 x comb_15
hk1(1:k1)=1/k1;
hk2=conv(hk1,hk1);
hk4=conv(hk2,hk2);
h1=conv(hk4,hk4);      
n11=ceil(length(h1)/d1);
n12=ceil(n/d1);

% Convolve sinusoidal input with impulse response of comb filter
y1=conv(x,h1);
y1=y1(1:n);
w1=kaiser(n,20);               % window type for fft 
ny1=norm(w1,2)*sqrt(n/2);
sy1=abs(fft(y1'.*w1))/ny1;
sy1db=20*log10(sy1);

% decimate impulse response and filter output by d1
for i=1:n11
  h1d(i)=d1*h1(d1*(i-1)+1);  
end
for i=1:n12
  y1d(i)=d1*y1(d1*(i-1)+1);  
end

h10=h1d;                  
len1=length(h10);
h10(len1+1:n)=0;	     % zero padded impulse response of comb filter
s1=abs(fft(h10));
s2db=20*log10(s1);

% compute response before decimation by d1

h1x=h1;                      % response before decimation by d1
len1=length(h1x);
h1x(len1+1:n)=0;	     % zero padded impulse response of comb filter
s1x=abs(fft(h1x));
s1db=20*log10(s1x);


% Compute equalizer coeffcients

fs3=fs/d1;
N3=2*m3-1;              % equalizer filter order
w3=pi*fp/fs3*z3;        % frequency of ith gain value
k3=round(n*fp/fs3*z3);  % index of ith gain value

for i=1:m3
  c3(i)=1/s1(k3(i));    % inverse passband gain values
end

% compute filter coefficients

for i=1:m3
  for j=1:m3
    B3(i,j)=2*cos((N3+2-2*j)*w3(i));
  end
end
h3=inv(B3)*c3';

for i=1:m3
  h3(m3+i)=h3(m3+1-i);
end

% round equalizer coefficients to b3 bits
h3=h3*64/max(abs(h3));   % equate maximum coeffcient to 64
h3=round(b0*h3)/b0;

%  Conversion to 2's complement
le3=ceil(length(h3)/2);
h31=h3(1:le3)*b0;
for i=1:le3
  if h31(i)<0
  h31(i)=128*b0+h31(i);
  end
end

% conversion to signed magnitude
h32=abs(h3(1:le3))*b0*2;

h3_bin=dec2bin(h32)
%h3_hex=dec2hex(h32)

% Convolve h1d with h3

h33=conv(h1d,h3);

h30=h33;
len3=length(h30);
h30(len3+1:n)=0;	% zero padded impulse response of equalized filter
s3=abs(fft(h30));
s3db=20*log10(s3);

%  convolve output of comb filter with equalizer impulse response

y2d=conv(y1d,h3);
y2d=y2d(1:n2)/d1;
w2=kaiser(n2,20);                       % window type for fft 
ny2=norm(w2,2)*sqrt(n2/2);
sy2=abs(fft(y2d'.*w2))/ny2;
sig_power_dB=10*log10(sum(sy2.^2)/2)    % signal power
sy2db=20*log10(sy2);

% plot filter response
figure(1)

subplot(2,1,1)
plot(d1*f1,sy1db)
axis([0 fs/d1 -180 0])
title('Comb filter Output before decimation')
grid

subplot(2,1,2)
plot(f2,sy2db)
axis([0 fs/(2*d1) -180 0])
title('Equalized Comb Filter Output')
grid

figure(2)               % Plot Impulse Response

subplot(2,1,1)
maxh1=max(h1);
impz(h1)
title('Impulse Response of Comb Filter Cascade')
axis([0 N1 0 1.1*maxh1])
grid;

subplot(2,1,2)
maxh3=max(h3);
minh3=min(h3);
impz(h3)
title('Impulse Response of Equalizer')
axis([-.01 2*m3-1 1.1*minh3-.1 1.1*maxh3])
grid;

figure(3)           % Plot Filter Response

subplot(3,1,1)
plot(d1*f1,s1db)
s1max=max(s1db);
axis([0 d1*fs3/2 -200 .5+s1max])
title('8-fold Cascade of Comb15 before Decimation by 10')
xlabel('Frequency [kHz]');
ylabel('Amplitude [dB]');
grid;

subplot(3,1,2)
plot(f1,s2db)
s1max=max(s1db);
axis([0 fs3/2 -60 .5+s1max])
title('8-fold Cascade of Comb15 after Decimation by 10')
xlabel('Frequency [kHz]');
ylabel('Amplitude [dB]');
grid;

subplot(3,1,3)
plot(f1,s3db)
s3max=max(s3db);
%axis([0 fs3/2 -1.2*s3max 1.2*s3max])
axis([0 fs3/2 -3 .4+s3max])
title(['Equalized Passband Response of 8-fold Cascade of Comb15'])
xlabel('Frequency [kHz]');
ylabel('Amplitude [dB]');
grid;

end

