[Hd]=peakingEQ(fc,g,Q)

function [Hd]=peakingEQ(fc,g,Q)

fs=44100;
% coefficient calculation
g=10^(g/20);
t0=2*pi*fc/fs;
if g >= 1
    beta=t0/(2*Q);
else
    beta=t0/(2*g*Q);
end    
a2=-0.5*(1-beta)/(1+beta);
a1=(0.5-a2)*cos(t0);
b0=(g-1)*(0.25+0.5*a2)+0.5;
b1=-a1;
b2=-(g-1)*(0.25+0.5*a2)-a2;

% SOS Form
b=2*[b0 b1 b2];
a=[1 -2*a1 -2*a2];
[sos,g]=tf2sos(b,a);
Hd=dfilt.df1sos(sos,g);