%%
% by Ljiljana Milic
% Supplemental material for Chapter X:
%% 10. Frequency Response Masquing Techniques
% Table of Contents
% 10. Frequency Response Masquing Techniques
%     10.1 NARROW BAND FILTERS
%         Example 10.1
%         Example 10.2
%     10.2 ARBITRARY BANDWIDTH DESIGN
%     10.2.1 Filter Synthesis Based on a FIR Complementary Filter Pair and
%      Two Masking Filters
%         Example 10.3
%     10.2.2 Filter Synthesis Based on a Complementary IIR Filter Pair and
%      Two FIR Masking Filters
%         Example 10.4
%%
%% 10.1 NARROW BAND FILTERS
%
% The narrowband filter H\left(z\right) is obtained as a cascade of periodic model filter G\left(z^M \right)and masking filter F\left(z\right).
%
% H\left(z\right)=G\left(z^M \right)F\left(z\right)
%% Example 10.1
% Specifications for lowpass filter:
% Passband cutoff: F_p =100\;Hz
% Stopband cutoff: F_s =150 Hz
% Sampling frequency: F_0 =2000 Hz
% Passband ripple \delta_p =0\ldotp 01, and the stopband ripple \delta_s =0\ldotp 01(a_s =40dB)
% Solution:
% For the frequency response masking approach we choose M = 4, thus we design the model filter G\left(z\right) with F_p =100\;Hz, and sampling frequency F_0 /M=500 Hz .
%
%  Design od the model filter G\left(z\right)

clear all, close all
[N,fk,mk,w]=firpmord([100,150],[1,0],[0.01,0.01],500);
fg = [0,100/250,150/250,1]; mk=[1,1,0,0];
g = firhalfband(22,100/250);
[G,w1]=freqz(g,1,1000);
figure (1)
plot(500*w1/(2*pi),20*log10(abs(G))), grid
xlabel('Frequency (Hz)'),ylabel('Gain(dB)'), title('Model filter')
axis([0,250,-60,5])

% Generation of the periodic model filter G\left(z^M \right)

gm = zeros(size(1:4*length(g)));
gm([1:4:length(gm)]) = g;
[Gm,w1]=freqz(gm,1,1000);
figure (2)
plot(2000*w1/(2*pi),20*log10(abs(Gm))),grid
xlabel('Frequency (Hz)'),ylabel('Gain(dB)')
axis([0,1000,-60,5])
title('Periodic model filter')

% Design of the masking filter F\left(z\right)

f0=[0,100/1000,375/1000,625/1000,875/1000,1]; m0=[1,1,0,0,0,0];
w0=[w',w(2)];
%ff=firpm(18,f0,m0,w0);
ff = firpm(18,f0,m0,[1,0.5,0.5]);
[F,w1]=freqz(ff,1,1000);
figure (3)
plot(2000*w1/(2*pi),20*log10(abs(F))), grid
xlabel('Frequency (Hz)'),ylabel('Gain(dB)'),title('Masking filter')
axis([0,1000,-60,5])

% Narrowband filter H\left(z\right)

H=F.*Gm;
figure (4)
plot(2000*w1/(2*pi),20*log10(abs(H))),grid
xlabel('Frequency (Hz)'),ylabel('Gain(dB)')
title('Narrow-band filter')
axis([0,1000,-60,5])
g = axes('Position',[0.35 0.65 0.50 0.1]);
plot(2000*w1(1:102)/(2*pi),abs(H(1:102))), axis([0,100,0.99,1.011])
disp('END OF EXAMPLE 10.1')

%
%% Example 10.2
% This example illustrates the application of a simple EMQF (Elliptic Minimal Q factors Filter) in constructing the lowpass narrowband IIR filter.
% We show a solution for the narrowband filter based on the identical EMQF subfilters where the same protype is used for the model and for the masking filter, i.e. F\left(z\right)=G\left(z\right), as illustrated in the Figure bellow.
%
%
% The narowband filter transfer functionH\left(z\right) is obtained as a product of the periodic model filter G\left(z^M \right)and the masking filter G\left(z\right):
% H\left(z\right)=G\left(z^M \right)G\left(z\right)
% where the transfer function G\left(z\right) is expressible as a sum of two stable allpass functions A_0^G \left(z\right) and A_1^G \left(z\right)
%
% G\left(z\right)=\frac{1}{2}\left(A_0^G \left(z\right)+A_1^G \left(z\right)\right)
% A_0^G \left(z\right)=\frac{\beta_1 -0\ldotp 5\left(1+\beta_1 \right)z^{-1} +z^{-2} }{1-0\ldotp 5\left(1+\beta_1 \right)z^{-1} +\beta_1 z^{-2} }
% A_1^G \left(z\right)\left(z\right)=\frac{\alpha_1 +z^{-1} }{1+\alpha_1 z^{-1} }\times \frac{\beta_2 -0\ldotp 5\left(1+\beta_2 \right)z^{-1} +z^{-2} }{1-0\ldotp 5\left(1+\beta_2 \right)z^{-1} +\beta_2 z^{-2} }
% Filter specifications
% Stopband cutoff: \omega_s =0\ldotp 1\pi
% Passband cutoff: \omega_p =0\ldotp 065\pi
% Minimal stopband attenuationa_s =30dB
% Solution:
% For G\left(z\right) we choose the 5th order  EMQF filter with the 3 dB cutt off at \omega_{3\textrm{dB}} =\left(1/3\right)\pi and the stopband edge frequency \omega_s =0\ldotp 4\pi. This choice guaranties  the stopband attenuation of 30 dB, and very small passband ripple with the passband cutoff at \omega_p =0\ldotp 2738\pi.
% The two allpass branches of the 5th order EMQF filter, A_0^G \left(z\right) and A_1^G \left(z\right), are implemented with only 5 multiplication constants \alpha_1, \beta_1, and \beta_2:
% A_0^G \left(z\right)=\frac{\beta_1 -0\ldotp 5\left(1+\beta_1 \right)z^{-1} +z^{-2} }{1-0\ldotp 5\left(1+\beta_1 \right)z^{-1} +\beta_1 z^{-2} }
% A_1^G \left(z\right)=\frac{\alpha_1 +z^{-1} }{1+\alpha_1 z^{-1} }\times \frac{\beta_2 -0\ldotp 5\left(1+\beta_2 \right)z^{-1} +z^{-2} }{1-0\ldotp 5\left(1+\beta_2 \right)z^{-1} +\beta_2 z^{-2} }

clear, close all

% Design of EMQF filter G\left(z\right)
% The function file emqf_LP.m is used for the EMQF design

N = 5; alpha = -0.5; Fs = 0.4; % Input parameters for emqf_LP.m
[A1b,A1i,A0b,A0i,beta,Ap,As,Fp,F3db]=emqf_igi(N,alpha,Fs/2); % Filter design
A=(conv(A1b,A0i)+conv(A0b,A1i))/2; % Numerator of G(z)
B=conv(A1i,A0i); % Denominator of G(z)
[G,f]=freqz(A,B,1024,2); 
figure (1)
plot(f,20*log10(abs(G))),ylabel('Gain, dB'),grid
xlabel('\omega/\pi')
title('EMQF filter, used as model and masking filter ')
axis([0,1,-50,2])

% Generation of the periodic model filter G\left(z^M \right),\;M=4

AA=zeros(size(1:4*length(A)-3));
AA(1:4:length(AA))=A;
BB=zeros(size(1:4*length(B)-3));
BB(1:4:length(BB))=B;
[Gm,f]=freqz(AA,BB,1024,2);
figure (2)
plot(f,20*log10(abs(Gm))),ylabel('Gain, dB'),grid
xlabel('\omega/\pi')
title('Periodic EMQF filter')
axis([0,1,-50,2])

% Computing frequency response of the narrowband filter H\left(z\right)

H=G.*Gm;
figure (3)
plot(f,20*log10(abs(H))),xlabel('\omega/\pi'),ylabel('Gain, dB'),grid
title('Narrowband filter')
axis([0,1,-50,2])
g(1)=axes('Position',[0.30 0.70 0.50 0.10]);
plot(f(1:74),abs(H(1:74))), axis([0,74/1024,0.9992,1.0001])
disp('Notice: H(z) needs 6 multiplications and 2 binary shifts.')
disp('END OF EXAMPLE 10.2')

%
%% 10.2 ARBITRARY BANDWIDTH DESIGN
%
%% 10.2.1 Filter Synthesis Based on a FIR Complementary Filter Pair and
%%  Two Masking Filters
%
% This method is based on the complementary pair of periodic filters and two masking filters as shown in the Figure bellow. The transfer function of the overall filter H\left(z\right) is represented in the form:
% H\left(z\right)=G\left(z^M \right)F_0 \left(z\right)+G_c \left(z^M \right)F_1 \left(z\right)
% where M is a positive integer,  G\left(z^M \right) is periodic model filter, and G_c \left(z^M \right)periodic complementary model filter, F_0 \left(z\right) and F_1 \left(z\right) are masking filters.
%
%
%%
%% Example 10.3
% This example demonstrates the filter synthesis method based on the complementary FIR filter pair, and the two masking filters according to the structure shown above.
% Specisications for the  wideband linear-phase FIR lowpass filter:
% Passband cutoff: F_p =600\;Hz
% Stopband cutoff: F_s =650 Hz
% Sampling frequency: F_0 =2000 Hz
% Maximal passband ripple a_p =\pm \;0\ldotp 1dB
% Minimal stopband attenuation a_s =40dB
%

clear all, close all

% Model filter design

Nord = 22; % Filter order
fp= 0.4; % Model filter cutoff
g = firhalfband(Nord,fp);              
[G,fr] = freqz(g,1,1024,2);

figure (1)
plot(fr,20*log10(abs(G))); grid
xlabel('Normalized frequency, \omega/\pi'),ylabel('Gain(dB)')
title('Model filter')
axis([0,1,-60,5])

% Design parameters for masking filters

ap = 0.1; as = 40;
delta_p = (10^(ap/20)-1)/(10^(ap/20)+1); % Passband ripple
delta_s = 10^(-as/20);                   % Stopband ripple

% Designing the masking filter F_0 \left(z\right)

freq = [0,0.6,0.85,1]; m0 = [1,1,0,0]; w = [delta_p/delta_s,1]; 
f0 = firpm(18,freq,m0,w);     
[F0,fr] = freqz(f0,1,1024,2);

%  Designing the masking filter F_1 \left(z\right)

freq = [0,0.4,0.625,1]; m1=[1,1,0,0]; 
f1 = firpm(20,freq,m1,w);     
[F1,fr] = freqz(f1,1,1024,2);

figure (2)
plot(fr,20*log10(abs(F0)),fr,20*log10(abs(F1)),'r'),grid 
xlabel('Normalized frequency, \omega/\pi'),ylabel('Gain(dB)')
title('Masking filters')
legend('Location','best')
legend('Filter F_0(z)','Filter F_1(z)')
axis([0,1,-60,5])

% Generation of the periodic model filter G\left(z^M \right)

gm = zeros(size(1:4*length(g)));        
gm([1:4:length(gm)]) = g;                  
[Gm,fr] = freqz(gm,1,1024,2);
figure (2)
plot(fr,20*log10(abs(Gm)))
hold on
xlabel('Normalized frequency, \omega/\pi'), ylabel('Gain(dB)')
title('Complementary periodic model filter pair')
axis([0,1,-60,5])

% Generation of the complementary periodic model filter

d = zeros(size(gm)); d(45) = 1;          % Delay branch
gmc = d - gm;                            % Complementary periodic model filter
[Gmc,fr]=freqz(gmc,1,1024,2);
plot(fr,20*log10(abs(Gmc)),'--b'), grid

% Implementation of the wide-band filter H\left(z\right)

f0=[0,f0,0]; % Equalization of the masking filters' lengths
h = conv(gm,f0) + conv(gmc,f1); % Impulse response of the overall filter
[H,Fr]=freqz(h,1,1024,2000);
figure (4)
plot(Fr,20*log10(abs(H))); 
xlabel('Frequency [Hz]'),ylabel('Gain, dB')
hold on
plot([0,650,650,1000],[0,0,-40,-40],'r')
title('Wide-band filter')
axis([0,1000,-60,5])
g(1)=axes('Position',[0.20 0.55 0.40 0.07]);
plot(Fr(1:650),20*log10(abs(H(1:650)))),axis([0,600,-0.1,0.1])
disp('END OF EXAMPLE 10.3')

%
%% 10.2.2 Filter Synthesis Based on a Complementary IIR Filter Pair and
%%  Two FIR Masking Filters
%
% The block diagram shown bellow represents the structure of the IIR filter synthsized using the frequency response masking technique.
%
%%
% The complementary model filter pair G\left(z\right)\left\lbrack G_c \left(z\right)\right\rbrack is expresible as a sum/difference of two stable allpass functionsA_0^G \left(z\right) and  A_1^G \left(z\right):
% G\left(z\right)\left\lbrack G_c \left(z\right)\right\rbrack =\frac{1}{2}\left(A_0^G \left(z\right)\pm A_1^G \left(z\right)\right)
% and consequently the overall filter transfer function H\left(z\right) is given by
%
% H\left(z\right)=\frac{A_0^G \left(z^M \right)+A_1^G \left(z^M \right)F_1 \left(z\right)}{2}F_0 \left(z\right)+\frac{A_0^G \left(z^M \right)-A_1^G \left(z^M \right)}{2}F_1 \left(z\right)
% where F_0 \left(z\right) and F_1 \left(z\right) are masking filters.
%
%% Example 10.4
% This example demonstrates the filter synthesis method based on the complementary IIR filter pair, and the two masking filters according to the structure shown above.
%
% Specisications for the  wideband IIR lowpass filter:
% Passband cutoff: F_p =600\;Hz
% Stopband cutoff: F_s =650 Hz
% Sampling frequency: F_0 =2000 Hz
% Maximal passband ripple a_p =\pm \;0\ldotp 1dB
% Minimal stopband attenuation a_s =40dB

clear all, close all

% Allpass branches

[b,a,z,p,k] = halfbandiir(5,0.4);  % Halfband filter design
beta2=(abs(p(2)))^2; beta3=(abs(p(4)))^2; % Coefficients of all-pass branches

% Complementary model filter pair G\left(z\right),G_c \left(z\right)

[A0,w1]=freqz([beta2,0,1],[1,0,beta2],512);  
[A1,w1]=freqz([0,beta3,0,1],[1,0,beta3],512);          
G=(A0+A1)/2; Gc=(A0-A1)/2; 

% Design of masking filters F_0 \left(z\right) and F_1 \left(z\right)

ap = 0.1; as = 40;
delta_p = (10^(ap/20)-1)/(10^(ap/20)+1); % Passband ripple
delta_s = 10^(-as/20);                   % Stopband ripple
freq = [0,0.6,0.85,1]; m0 = [1,1,0,0]; w = [delta_p/delta_s,1]; % Seting design parameters for masking filter f0
f0 = firpm(18,freq,m0,w);     % Designing the masking filter f0
f0 = [0,f0,0];
freq = [0,0.4,0.625,1]; m1=[1,1,0,0];  % Seting design parameters for masking filter f1
f1 = firpm(20,freq,m1,w);     % Designing the masking filter f1
[F0,w1]=freqz(f0,1,512);
[F1,w1]=freqz(f1,1,512);

% Generation of periodic model filter pair G\left(z^M \right),G_c \left(z^M \right)

[A0m,w1]=freqz(upsample([beta2,0,1],4),upsample([1,0,beta2],4),512);  
[A1m,w1]=freqz(upsample([0,beta3,0,1],4),upsample([1,0,beta3],4),512); 
Gm=(A0m+A1m)/2;    Gcm=(A0m-A1m)/2;          

% Implementation of the wideband filter H\left(z\right)

H=Gm.*F0+Gcm.*F1;
figure (1)
plot(2000*w1/(2*pi),20*log10(abs(H))); 
xlabel('Frequency [Hz]'),ylabel('Gain, dB')
title('Wideband filter')
axis([0,1000,-60,5])
hold on
plot([0,650,650,1000],[0,0,-40,-40],'r')
Fr=2000*w1/(2*pi);
g(1)=axes('Position',[0.20 0.55 0.40 0.07]);
plot(Fr(1:325),20*log10(abs(H(1:325)))),axis([0,600,-0.1,0.1])
disp('END OF EXAMPLE 10.4')
disp('                      END OF CHAPTER X')

%

function [A1b,A1i,A0b,A0i,beta,Ap,As,Fp,F3db]=emqf_igi(N,alpha,Fs)
%
% emqf_0	EMQF (Elliptic Minimal Q Factors) digital filter design.
%		          [A1b,A1i,A0b,A0i,beta,Ap,As,Fp,F3db]=emqf_LP(N,alpha,Fs) designs an N'th order
%               lowpass digital EMQF filter with normalized stopband edge
%               frequency Fs. The sampling frequency is assumed to be F0=1.
%               The constant alfa is chosen by designer. For given Fp and Fs
%               (Fp and Fs the boundary frequencies of the specificatin),
%               the permited values alfa are in the range:
%
%   -(1-(tan(pi*Fa)^2)/(1+(tan(pi*Fa)^2)<alfa<-(1-(tan(pi*Fp)^2)/(1+(tan(pi*Fp)^2).
%
%	   The value of alfa approximates 
%              
%               (1-(tan(pi*Fa)(tan(pi*Fp))/(1+(tan(pi*Fa)(tan(pi*Fp))
% 
%               and is adjusted to be represented as a simple combination
%               of powers of two. 
%    
%     In the case of a half-band filter alfa=0         
%
%       emqf_1 returns:
%              * Allpass branches in vectors A1b,A1i,A0b,A0i
%              * Normalized frequencies (F0=1) in vector f, length 500. 
%              * Square pole magnitudes of the complex poles in vector beta. 
%              * passband ripple Ap (dB)
%              * minimal stopband attenuation As (dB)
%              * passband boundary frequency of EMQF filter Fp
%              * 3 dB attenuation frequency F3db
%
%

disp(' EMQF design')
disp('===============================')
disp([sprintf(' N     = %6u', N)])
disp([sprintf(' alpha = %15.9f', alpha)])
disp([sprintf(' F_s   = %15.9f',2*Fs)])

F3db = (1/pi)*atan(sqrt((1+alpha)/(1-alpha)));

neven = N-2*fix(N/2);
if neven < 0.5
   disp('------- error -----------------')
   disp('N should be odd number')
   disp('choose new N = 3, 5, 7, 9, 11, ...')
   disp('-------------------------------')
   H=0; w=0; zzd=0; ppd=0; Ap=0; As=0; Fp=0; F3db=0;omega=0;co=0;
   return
end

if F3db >= Fs
   disp([sprintf(' F_3dB = %15.9f', F3db)])
   disp('------- error -----------------')
   disp('Fs should be larger than F3dB')
   disp(['choose new Fs    > ' num2str(F3db) '  or'])
   disp(['choose new alpha < ' num2str(-cos(2*pi*Fs))])
   disp('-------------------------------')
   H=0; w=0; zzd=0; ppd=0; Ap=0; As=0;omega=0;co=0;
   return
end

Fp = (1/pi)*atan(((1+alpha)/(1-alpha))*(1/tan(pi*Fs)));
disp([sprintf(' F_p   = %15.9f', 2*Fp)])
disp([sprintf(' F_3dB = %15.9f', 2*F3db)])

xi = tan(pi*Fs)/tan(pi*Fp);
if xi >= sqrt(2)
   t = .5*(1-(1-1/xi^2)^(1/4))/(1+(1-1/xi^2)^(1/4));
   q = t+2*t^5+15*t^9+150*t^13;
else
   t = .5*(1-1/sqrt(xi))/(1+1/sqrt(xi));
   qp = t+2*t^5+15*t^9+150*t^13;
   q = exp(pi^2/log(qp));
end
L  = (1/4)*sqrt(1/q^N-1);
Ap = 10*log10(1+1/L);
As = 10*log10(1+L);
disp([sprintf(' A_p   =     %15.9g', Ap)])
disp([sprintf(' A_s   = %15.9f', As)])

a = xi;
b = 1; 
for ind = 1:(N-1)/2
  x = ellipj( ((2*ind-1)/N + 1)*ellipke(1/a^2), 1/a^2);
  kor = 2*sqrt(a)*tan(pi*Fp)/(1+a*tan(pi*Fp)^2);
  b(ind) = ((a+x^2) - kor*sqrt((1-x^2)*(a^2-x^2)) )/...
           ((a+x^2) + kor*sqrt((1-x^2)*(a^2-x^2)) );
end
a1 = -(1-tan(pi*F3db))/(1+tan(pi*F3db));
beta = sort(b);
disp([sprintf(' alpha1 = %15.9f', a1)])
disp(' beta = ')
disp(beta')

p1 = [1 a1];
for ind = 2:2:(N-1)/2
  p1 = conv(p1,[1 alpha*(1+beta(ind)) beta(ind)]);
end
p2 = 1;
for ind = 1:2:(N-1)/2
  p2 = conv(p2,[1 alpha*(1+beta(ind)) beta(ind)]);
end
%l =1
A1b=p1;
A1i=fliplr(A1b);
A0b=p2;
A0i=fliplr(A0b);
end