% a simple demo for parallel filter design
%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
[Bm,Am,FIR]=parfiltdesfr(w(start:stop),spH(start:stop),p,NFIR);
filtH=parfiltfresp(Bm,Am,FIR,w); % freq. response of the filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plotting
button=1;
while button~=2,
semilogx(fr,db(spH),'b','LineWidth',2);
hold on;
semilogx(fr,db(filtH),'r','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),':');
if NFIR>0, % compute and plot the response of the FIR part
HFIR=freqz(FIR,1,w);
semilogx(fr,db(HFIR),'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,20,'Blue: target, Red: filter, Dotted: sections','FontName','Times','Fontsize',SIZE);
s=sprintf('Total filter order: %i', length(p));
text(20,25,s,'FontName','Times','Fontsize',SIZE);
axis([18 22000 -30 30]);
[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
[Bm,Am,FIR]=parfiltdesfr(w(start:stop),spH(start:stop),p,NFIR);
filtH=parfiltfresp(Bm,Am,FIR,w); % freq. response of the filter
end;