[Hd]=highpassEQ(fc,g,N)

function [Hd]=highpassEQ(fc,g,N)

fs=44100;
% parameter limits 
if fc < 20
    fc=20;
end
if fc > 20000
    fc=20000;
end

% coefficient calculation
if N==1
    g=10^(g/20);
    a1=-0.5*tan(pi*(fc/fs-0.25));
    b0=g*(0.25+0.5*a1);
    b1=-b0;
    b2=0;
    a2=0;
    % SOS Form
    b=2*[b0 b1 b2];
    a=[1 -2*a1 -2*a2];
    [sos,g]=tf2sos(b,a);
    Hd=dfilt.df1sos(sos,g);
elseif N==2
    g=10^(g/20);
    t0=tan(pi*fc/fs);
    t1=1+sqrt(2)*t0+t0^2;
    a1=(1-t0^2)/t1;
    a2=0.5*(sqrt(2)*t0-1-t0^2)/t1;
    b0=g/(2*t1);
    b1=-2*b0;
    b2=b0;
    % SOS Form
    b=2*[b0 b1 b2];
    a=[1 -2*a1 -2*a2];
    [sos,g]=tf2sos(b,a);
    Hd=dfilt.df1sos(sos,g);
elseif N==3
    g=10^(g/20);
    t0=tan(pi*fc/fs);
    t1=1+t0+t0^2;
    a01=0.5*(1-t0)/(1+t0);
    a11=(1-t0^2)/t1;
    a12=(t0-t0^2-1)/(2*t1);
    b00=g/2*1/(t0+1);
    b01=-b00;
    b10=1/(2*t1);
    b11=-2*b10;
    b12=b10;
    a02=0;
    b02=0;
    % SOS Form
    b0=2*[b00 b01 b02];
    a0=[1 -2*a01 -2*a02];
    [sos0,g0]=tf2sos(b0,a0);
    Hd0=dfilt.df1sos(sos0,g0);
    b1=2*[b10 b11 b12];
    a1=[1 -2*a11 -2*a12];
    [sos1,g1]=tf2sos(b1,a1);
    Hd1=dfilt.df1sos(sos1,g1);
    Hd=dfilt.cascade(Hd0,Hd1);
elseif N==4
    g=10^(g/20);
    t0=tan(pi*fc/fs);
    t1=1+2*cos(1/8*pi)*t0+t0^2;
    t2=1+2*cos(3/8*pi)*t0+t0^2;    
    a01=(1-t0^2)/t1;
    a02=(2*cos(1/8*pi)*t0-1-t0^2)/(2*t1);
    a11=(1-t0^2)/t2;
    a12=(2*cos(3/8*pi)*t0-1-t0^2)/(2*t2);
    b00=0.5*1/t1;
    b01=-2*b00;
    b02=b00;
    b10=g/(2*t2);
    b11=-2*b10;
    b12=b10;
    % SOS Form
    b0=2*[b00 b01 b02];
    a0=[1 -2*a01 -2*a02];
    [sos0,g0]=tf2sos(b0,a0);
    Hd0=dfilt.df1sos(sos0,g0);
    b1=2*[b10 b11 b12];
    a1=[1 -2*a11 -2*a12];
    [sos1,g1]=tf2sos(b1,a1);
    Hd1=dfilt.df1sos(sos1,g1);
    Hd=dfilt.cascade(Hd0,Hd1);
end