%**************************************************************************
% Main parameters
%**************************************************************************
N_pred = 500;

%**************************************************************************
% Auxillary parameters
%**************************************************************************
N_FFT     =   512*4;
f_s       = 44100;
LineWidth =     1;

%**************************************************************************
% Desired frequency response
%**************************************************************************
H_des                   = zeros(N_FFT/2+1,1)- 30; 
H_des(1)                = 30;
H_des(500:600)          = 40;
beta                    = 0.95;
H_des                   = filtfilt(1-beta,[1 -beta],H_des);
H_des_mag               = 10.^(H_des/20);

%**************************************************************************
% Minimum-phase FIR filter
%**************************************************************************
S_hh                     = 1./[H_des_mag; conj(H_des_mag(end-1:-1:2))].^2;
s_hh                     = ifft(S_hh);
[b_fir_minphas, pow_err] = levinson(s_hh,N_pred-1);
b_fir_minphas            = b_fir_minphas / sqrt(pow_err);
a_fir_minphas            = 1;

x = randn(4000,1)*2;
y = filter(b_fir_minphas,a_fir_minphas,x);

%**************************************************************************
% Result analysis
%**************************************************************************
S_xx = pwelch(x,hanning(N_FFT),N_FFT*(1-1/16),N_FFT);
S_yy = pwelch(y,hanning(N_FFT),N_FFT*(1-1/16),N_FFT);

f = 0:1/N_FFT:0.5;
f = f * 2;

fig = figure(2);
set(fig,'Units','Normalized');
set(fig,'Position',[0.1 0.1 0.8 0.5]);

subplot('Position',[0.1 0.12 0.8 0.4]);
plot(f,H_des,'r', ...
     f,10*log10(S_yy),'b','LineWidth',LineWidth);
axis([0 1 -45 45])
xlabel('Normalized frequency \Omega/\pi')
grid on;

ylabel('dB')
legend('Desired magnitude response', ...
       'Estimated power speactral dessity of the filter output'); 
   
subplot('Position',[0.1 0.57 0.8 0.4]);
plot(f,H_des,'r', ...
     f,10*log10(S_xx),'b','LineWidth',LineWidth);
axis([0 1 -45 45])
grid on;
set(gca,'XTickLabel','');
ylabel('dB')
legend('Desired power spectral desnity', ...
       'Estimated power speactral dessity of the filter input');

  