% ----------------------------------------------------------------------------
% Example output from hpeq.m
%
% Elliptic bandstop filter
% ----------------------------------------------------------------------------
% Sophocles J. Orfanidis
% ECE Department, Rutgers University
% 94 Brett Road, Piscataway, NJ 08854-8058
%
% Email:   orfanidi@ece.rutgers.edu
% Date:    June 15, 2005
%
% Reference: Sophocles J. Orfanidis, "High-Order Digital Parametric Equalizer
%            Design," J. Audio Eng. Soc., vol.53, pp. 1026-1046, November 2005.
%
% Web Page: http://www.ece.rutgers.edu/~orfanidi/hpeq
%
% tested with MATLAB R11.1 and R14
% ----------------------------------------------------------------------------

clear all;

fs=40; f1=12; f2=14; f1s=11.5; f2s=14.5;
G0 = 0; G = -Inf; GB = -30; type = 3; Gs = -0.1; tol=eps;

w1=2*pi*f1/fs; w2=2*pi*f2/fs; w1s=2*pi*f1s/fs; w2s=2*pi*f2s/fs;

[w0,Dw] = bandedge(w1,w2,1); Dws = w2s-w1s;

Nexact = hpeqord(G0, G, GB, Gs, Dw, Dws, type),
N=ceil(Nexact)

[B,A,Bh,Ah] = hpeq(N, G0, G, GB, w0, Dw, type, Gs, tol)

f = linspace(0,20,1001); w = 2*pi*f/fs;

H = abs(fresp(B,A,w)).^2;

% H = abs(fresp(Bh,Ah,w,w0));  % alternative computation of frequency response

Gb = -3; Dwb = bandwidth(N,G0,G,GB,Gb,Dw,type,Gs);                         % 3-dB width
[w1b,w2b] = bandedge(w0,Dwb); f1b = fs * w1b/2/pi; f2b = fs * w2b/2/pi;    % 3-dB frequencies

figure;
plot(f,H,'r-', [f1,f2,f1b,f2b,f1s,f2s], 10.^([GB,GB,Gb,Gb,Gs,Gs]/10), 'b.');
ylim([0,1.05]); ytick(0:0.25:1);
xlim([0,20]); xtick(0:2:20);
title('Elliptic, Bandstop');
xlabel('{\it f}  (kHz)'); ylabel('magnitude squared')
grid;
Nexact =

    4.0308


N =

     5


B =

    0.1572    0.1446    0.1572         0         0
    1.3739    2.4677    3.7550    2.4677    1.3739
    2.9724    5.4008    8.2874    5.4008    2.9724


A =

    1.0000    0.7325    0.5936         0         0
    1.0000    1.7559    2.6161    1.7075    0.9473
    1.0000    1.6466    2.2804    1.3911    0.7221


Bh =

    0.1572   -0.1572         0
    1.3739   -2.6207    1.3739
    2.9724   -5.8050    2.9724


Ah =

    1.0000   -0.5936         0
    1.0000   -1.8202    0.9473
    1.0000   -1.5823    0.7221