% ----------------------------------------------------------------------------
% Example output from hpeq.m
%
% Chebyshev type-1 & type-2, bandpass/bandstop filters
% ----------------------------------------------------------------------------
% 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;

disp('% ---------- Chebyshev type-1 --------------------------------------------------');

fs=40; f1=4; f2=6; f1s=3.5; f2s=6.5;
G0 = -Inf; G = 0; GB = -0.1; type = 1; Gs = -30;

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)

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);                            % 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('Chebyshev-1, Bandpass');
xlabel('{\it f}  (kHz)'); ylabel('magnitude squared')
grid;

disp('% ---------- Chebyshev type-2 --------------------------------------------------');

fs=40; f1=11.5; f2=14.5; f1s=12; f2s=14; w1=2*pi*f1/fs;
G0 = 0; G = -Inf; GB = -0.1; type = 2; Gs = -30;

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)

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);                            % 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('Chebyshev-2, Bandstop');
xlabel('{\it f}  (kHz)'); ylabel('magnitude squared')
grid;
% ---------- Chebyshev type-1 --------------------------------------------------

Nexact =

    6.1718


N =

     7


B =

    0.1079         0   -0.1079         0         0
    0.0124         0   -0.0248         0    0.0124
    0.0120         0   -0.0239         0    0.0120
    0.0117         0   -0.0234         0    0.0117


A =

    1.0000   -1.3512    0.8874         0         0
    1.0000   -2.7531    3.7908   -2.6810    0.9496
    1.0000   -2.7167    3.6699   -2.5218    0.8639
    1.0000   -2.7044    3.6149   -2.4285    0.8073


Bh =

    0.1079    0.1079         0
    0.0124    0.0248    0.0124
    0.0120    0.0239    0.0120
    0.0117    0.0234    0.0117


Ah =

    1.0000   -0.8874         0
    1.0000   -1.8456    0.9496
    1.0000   -1.7947    0.8639
    1.0000   -1.7776    0.8073

% ---------- Chebyshev type-2 --------------------------------------------------

Nexact =

    6.1718


N =

     7


B =

    1.1706    1.0931    1.1706         0         0
    0.4381    0.8181    1.2581    0.8181    0.4381
    0.5645    1.0543    1.6213    1.0543    0.5645
    0.9466    1.7678    2.7186    1.7678    0.9466


A =

    1.0000    0.5704    0.2216         0         0
    1.0000    1.7456    2.5389    1.6824    0.9324
    1.0000    1.6370    2.2172    1.4089    0.7558
    1.0000    1.3958    1.5124    0.8432    0.4082


Bh =

    1.1706   -1.1706         0
    0.4381   -0.8761    0.4381
    0.5645   -1.1291    0.5645
    0.9466   -1.8932    0.9466


Ah =

    1.0000   -0.2216         0
    1.0000   -1.7387    0.9324
    1.0000   -1.5062    0.7558
    1.0000   -0.9897    0.4082