% a simple demo for a parallel equalizer
% - similar to parfiltdemo, but this one equalizes the response
% instead of modeling it
%
% it starts with a few poles and then USE YOUR MOUSE!
% LEFT CLICK: put a pole at a certain frequency
% RIGHT CLICK: remove the pole closest to the click
% MIDDLE CLICK: quits the programme
%
% So just click with left where you want better resolution (where the fit
% is bad), and with right where you think you don't need it (the fit is too good)
%
% Please note that this is just for fun (and getting an idea on how the
% parallel filter works).
% Typically you can get more accurate response by the automatic pole 
% positioning methods. For the overview and comparison of various techniques, see
% 
% Balázs Bank, "Loudspeaker and room response equalization using parallel filters:
% comparison of pole positioning strategies," Proceedings of the 51st AES Conference,
% Helsinki, Finland, August 2013. 
%
% Matlab code and publications about parallel filters:
% http://www.mit.bme.hu/~bank/parfilt
%
% C. Balazs Bank, 2013.


% parameters:
PSTART=10; % number of pole frequencies to start with
SMOOTH=6; % fractional octave smoothing, e.g., SMOOTH=6 is 6th octave - 0 is no smoothing
PLOTSECT=1; % 1: plot the response of the separate sections, 0: don't plot


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computing smoothed loudspeaker-room response

Fs=44100;
load roomresp;

[cp,impresp]=rceps(5*impresp); %making mimumumphase
if SMOOTH>0,
    [fr,spH]=tfplots(impresp,'b',Fs,SMOOTH,'complex'); % smoothed loudspeaker-room response
else
    [fr,spH]=tfplot(impresp,'b',Fs); % un-smoothed loudspeaker-room response
end;
w=2*pi*fr/Fs;
% we will use the frequency range fr(start:stop) in the filter design
[tmp,start]=min(abs(fr-20)); % start frequency index
[tmp,stop]=min(abs(fr-Fs/2)); % stop frequency index

% a simple target respone
[Bh,Ah]=butter(4,30/(Fs/2),'high'); % target is a highpass
spT=freqz(Bh,Ah,w);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% logarithmic pole positioning as a first guess

fplog=logspace(log10(30),log10(20000),PSTART);
fp=fplog;
p=freqpoles(fp);
NFIR=0; % no FIR part is needed for a mimimum-phase filter

%parallel filter design - no FIR part in this case
[Bm,Am,FIR]=parfiltidfr(w(start:stop),spH(start:stop),spT(start:stop),p,NFIR);
filtH=parfiltfresp(Bm,Am,FIR,w); % freq. response of the filter
eqdH=filtH.*spH.';




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plotting


button=1;
while button~=2,

    


    semilogx(fr,db(spH)+0,'k--'); % original respolse
    hold on;


    semilogx(fr,db(spT),'b','LineWidth',2); % target
    semilogx(fr,db(eqdH),'r','LineWidth',2);
    
    semilogx(fr,db(filtH)-30,'b','LineWidth',2);

    if PLOTSECT, % compute and plot the response of the sections
        Z=exp(-j*w(:)); %creating z^-1
        Z2=Z.^2; %creating z^-2
        Hsectv=zeros(length(w),size(Am,2));
        for k=1:size(Am,2), %second-order sections
            Hsectv(:,k) = (Bm(1,k)+Bm(2,k)*Z)./(Am(1,k)+Am(2,k)*Z+Am(3,k)*Z2);
        end;
        semilogx(fr,db(Hsectv)-30,':'); 
        if NFIR>0, % compute and plot the response of the FIR part, if there is
            HFIR=freqz(FIR,1,w);
            semilogx(fr,db(HFIR)-30,'k');
        end;
    end;


    plot(fp,0,'ko','LineWidth',2); % plot pole frequencies
    hold off;
    


    SIZE=10;
    set(gca,'FontName','Times','Fontsize',SIZE);
    xlabel('Frequency [Hz]');
    ylabel('Magnitude [dB]');
    title('USE YOUR MOUSE: Left clik: add pole - Right click: remove - Middle click: finnish');
    text(500,10,'Black: original, Blue: target, Red: equalized','FontName','Times','Fontsize',SIZE);
    text(500,-20,'Blue: equalizer response, dotted: sections','FontName','Times','Fontsize',SIZE);
    s=sprintf('Total filter order: %i', length(p));
    text(20,15,s,'FontName','Times','Fontsize',SIZE);

    
    axis([18 22000 -60 20]);

    [x_cord,y_cord, button] = ginput(1); % ask for mouse click
    if button==1, % % if left click, add new pole at click
        fp=[fp x_cord];
        fp=sort(fp);
    end;
    if button==3, % if right click, delete the nearest pole to the click
        [val,ind]=min(abs(log(fp)-log(x_cord)));
        fp=[fp(1:ind-1) fp(ind+1:end)];
    end;
    p=freqpoles(fp);

    %parallel filter design - no FIR part in this case
    [Bm,Am,FIR]=parfiltidfr(w(start:stop),spH(start:stop),spT(start:stop),p,NFIR);
    filtH=parfiltfresp(Bm,Am,FIR,w); % freq. response of the filter
    eqdH=filtH.*spH.';


end;