Untitled

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

is obtained as a cascade of periodic model filter

and masking filter

Example 10.1

Specifications for lowpass filter:

Passband cutoff:

Stopband cutoff:

Sampling frequency:

Passband ripple

, and the stopband ripple

(

dB)

Solution:

For the frequency response masking approach we choose M = 4, thus we design the model filter

with

Hz, and sampling frequency

Hz .

Design od the model filter

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

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

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=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')

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.

, as illustrated in the Figure bellow.

The narowband filter transfer function

is obtained as a product of the periodic model filter

and the masking filter

where the transfer function

is expressible as a sum of two stable allpass functions

and

Filter specifications

Stopband cutoff:

Passband cutoff:

Minimal stopband attenuation

Solution:

For

we choose the 5th order EMQF filter with the 3 dB cutt off at

and the stopband edge frequency

. This choice guaranties the stopband attenuation of 30 dB, and very small passband ripple with the passband cutoff at

The two allpass branches of the 5th order EMQF filter,

and

, are implemented with only 5 multiplication constants

, and

clear, close all

Design of EMQF filter

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

 EMQF design
===============================
 N     =      5
 alpha =    -0.500000000
 F_s   =     0.400000000
 F_p   =     0.273837466
 F_3dB =     0.333333333
 A_p   =        0.0031204598
 A_s   =    31.437217336
 alpha1 =    -0.267949192
 beta = 
    0.3445
    0.7871

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

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=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.')

Notice: H(z) needs 6 multiplications and 2 binary shifts.

disp('END OF EXAMPLE 10.2')

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

is represented in the form:

where M is a positive integer,

is periodic model filter, and

periodic complementary model filter,

and

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:

Stopband cutoff:

Sampling frequency:

Maximal passband ripple

Minimal stopband attenuation

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

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

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

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

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')

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

is expresible as a sum/difference of two stable allpass functions

and

and consequently the overall filter transfer function

is given by

where

and

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:

Stopband cutoff:

Sampling frequency:

Maximal passband ripple

Minimal stopband attenuation

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

[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

and

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

[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=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')

END OF EXAMPLE 10.4

disp(' END OF CHAPTER X')

                      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