% Parallel second-order design example using a guitar admittance response
%
%   Details about the design algorithm can be found in the paper
%   Balazs Bank, "Direct Design of Parallel Second-order Filters for
%   Instrument Body Modeling", International Computer Music Conference,
%   Copenhagen, Denmark, Aug. 2007.
%
%   http://www.acoustics.hut.fi/go/icmc07-parfilt
%   http://www.mit.bme.hu/~bank/parfilt
%
%   C. Balazs Bank, Helsinki University of Technology, June 2007.


disp('Parallel second-order design example using a guitar admittance response');
disp('See comments for details, and the help files of the functions used.');

load guitadm;
Fs=44100;
[ceps,impresp]=rceps(impresp); % making it minimumphase

%%%%%%%%%%%%%%%%%%% THESE ARE THE MAIN PARAMETERS INFLUENCING THE DESIGN
lambda=0.65; %warping parameter
ORDER=50; %the order of the final filter, the number of second-order sections is ORDER/2
NFIR=1; %number of FIR coefficients - 1 means a simple feedforward path
%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Determining the pole set by a warped design (may take a while)...');
%determining the pole set by a warped design using the first 10000 samples
p=wppoles(impresp(1:10000),ORDER,lambda);
disp('Pole set ready.');

disp('Designing the parallel filter (this is faster)...');
%designing the parallel filter using the first 10000 samples of the impulse response
[Bm,Am,FIR]=parfiltdes(impresp(1:10000),p,NFIR);
disp('Filter ready.');

disp('Computing the impulse response of the designed filter...');
%calculating the impulse response
imp=zeros(10000,1);
imp(1)=1;
modresp=parfilt(Bm,Am,FIR,imp);
disp('Impulse response ready.');


%%%%%%%%%%%%%%%%%%%%%
%Plotting the results

figure(1); %impulse response
plot(impresp(1:10000),'k');
hold on;
plot(modresp,'r');
hold off;
axis([1 1000 -0.02 0.08]);
xlabel('Time [samples]');
ylabel('Amplitude');
title('Time-domain response - Black: original, Red: modeled');

figure(2); %logarithmically smoothed spectrum
tfplot(impresp(1:10000),'k'); %plotting the original response in logarithmic scale
hold on;
tfplot(modresp,'r'); %plotting the result in logaritmic scale
axis([20 20000 -30 30]);
hold off;
title('Frequency-domain response - Black: original, Red: modeled');
xlabel('Frequency [Hz]');
ylabel('Magnitude [dB]');