function filter=filterdesign(freq,ZQ,PQ,g);
% returns TF of filter
% frequencies are in f (not omega) (if freq='zpk', just returns zpk function)
% ZQ: matrix of zeros, Qs (put Q=0 for 1st order. otherwise, 2nd order)
% PQ: matrix of zeros, Qs (put Q=0 for 1st order. otherwise, 2nd order)
% g: filter gain
% Autor: Yuta Michimura

if length(ZQ)==0;
    zeros=[];
    Qzeros=[];
else
    zeros=ZQ(:,1)*2*pi;
    Qzeros=ZQ(:,2);
end
if length(PQ)==0;
    poles=[];
    Qpoles=[];
else
    poles=PQ(:,1)*2*pi;
    Qpoles=PQ(:,2);
end
z=[];
p=[];
for kk=1:length(Qzeros);
    if Qzeros(kk) ~= 0;
        a=1/(2*Qzeros(kk));
        b=sqrt(a^2-1);
        z=[z,-zeros(kk)*[a+b,a-b]];
    else
        z=[z,-zeros(kk)];
    end
end
for kk=1:length(Qpoles);
    if Qpoles(kk) ~= 0;
        a=1/(2*Qpoles(kk));
        b=sqrt(a^2-1);
        p=[p,-poles(kk)*[a+b,a-b]];
    else
        p=[p,-poles(kk)];
    end
end

F=zpk(z,p,g);
if freq(1) ~= 'z';
    [Fgain,Fphase,omega]=bode(F,2*pi*freq);
    filter=squeeze(Fgain.*exp(i.*Fphase/180*pi));
else
    filter=F;
end