function [p]=freqpoles(fr,Fs,Q);
% FREQPOLES - Pole set generation for parallel filter design.
% [P]=freqpoles(FR,Fs) creates poles from the given frequency vector
% FR with sampling frequency Fs such that the frequency responses
% of the second-order basis functions cross approximatelly at their
% -3 dB point.
% (See the equations in my 128th AES Convention paper.)
%
% [P]=freqpoles(FR,Fs,Q) creates poles from the given frequency vector
% FR with sampling frequency Fs but it uses the quality factors
% given in Q vector to set the radius of the poles. (-3dB bandwidth
% is DeltaF=FR/Q.)
%
% http://www.mit.bme.hu/~bank/parfilt
%
% C. Balazs Bank, 2010.
if nargin<2, %default sampling rate
Fs=44100;
end;
fr=fr(:); % making a column vector
wp=2*pi*fr/Fs; %discrete pole frequencies
if nargin<3, % we compute the bandwidth from the pole frequencies
wp=sort(wp); % sorting the frequencies
pnum=length(wp);
dwp=zeros(pnum,1);
for k=2:pnum-1,
dwp(k)=(wp(k+1)-wp(k-1))/2;
end;
dwp(1)=(wp(2)-wp(1));
dwp(pnum)=(wp(pnum)-wp(pnum-1));
end;
if nargin==3, % since Q are given, we compute the bandwidth from Q
Q=Q(:); % making a column vector
dwp=wp./Q;
end;
p=exp(-dwp/2).*exp(j*wp); % computing poles from center frequency wp and bandwidth dwp
p=[p; conj(p)]; %pole pairs one after the other
p=p(:); %making it a column vector