%%
% by Ljiljana Milic
% Supplemental material for Chapter VII:
%% 7. Lth-Band Digital Filters
%
% Table of Contents
% 7. Lth-Band Digital Filters
%     Introduction
%     7.1 Lth-Band  Linear-Phase FIR Filters
%         Example 7.1
%     7.2 Linear-Phase FIR Halfband Filters
%         Example 7.2
%         Example 7.3
%         Example 7.4
%     7.3 IIR Halfband Filter Design in MATLAB
%         Example 7.5
%     7.4 IIR Halfband Filters with Approximately Linear Phase
%         Example 7.6
%% Introduction
%
% The common characteristic of Lth-band lowpass filters is that the 6 dB (or 3dB) cutoff angular frequency is located at \pi /L, and the transition band is approximately symmetric around this frequency. In time domain, the impulse response of an Lth-band digital filter has zero valued samples at the multiples of L samples counted away from the central sample to the right and left directions. Actually, an Lth-band filter has the zero crossings at the regular distance of L samples thus satisfying the intersymbol interference property. The Lth-band filters are also called the Niquist filters.
%
%% 7.1 Lth-Band  Linear-Phase FIR Filters
%% Example 7.1
% This example presents a factor-of-3 interpolator constructed with the third-band filter (L = 3) of the length N = 21. The m-file firnyquist is used to compute impulse response coefficients. The example is aimed to demonstrate the following property of an Lth-band filter:  the sample values of the original signal are preserved in the interpolated signal, and in-between samples are the interpolated values.

close all, clear all

% 3rd-band filter design

h = firnyquist(20,3,0.3);
[H,F] = freqz(h,1,1000,2);

% Displayng the 3rd-band filter impulse response

figure 
stem((0:20),h)
xlabel('Time index n')
ylabel('Amplitude')
title('Impulse response')
axis([0,20,-0.1,0.4])

% Displaying the 3rd-band filter gain response

figure
plot(F,20*log10(abs(H))), grid
axis([0,1,-60,5])
xlabel('\omega/\pi')
ylabel('Gain, dB')
title('Gain Response')

% Generating the test signal x

n=0:16;
x=0.6*cos(2*pi*1.2*n/16)+0.8*sin(2*pi*0.7*n/16);

% Displaying the test signal

figure 
stem((0:16),x)
xlabel('Time index n')
ylabel('Amplitude')
title('Test signal')
axis([0,16,-1,1])

% Upsampling and interpolating the test signal

xu=upsample(x,3); % Upsampling-by-3
y=3*filter(h,1,xu); % Interpolating
y= 3*conv(h,xu);


% Displaying interpolated signal y

figure 
stem((10:3:length(xu)+7),x,'r')
hold on
stem(0:length(y)-1,y,'-.*k')
xlabel('Time index, n'), ylabel('y[m]')
legend('Exact values','Interpolated values')
axis([0,60,-1,1]);
title('Interpolated signal')
disp('In the interpolated signal y[m], every 3rd amplitude preserves') 
   disp('the value of the original test signal x[n], see red lines.')
disp('                    ')
disp('END OF EXAMPLE 7.1')

%
%% 7.2 Linear-Phase FIR Halfband Filters
%% Example 7.2
% In this example, we design a minimum order equaripple halfband FIR filter to satisfy the following specifications:
% The minimum stopband attenuation a_s =40\;\textrm{dB}(maximum stopband gain, \left.-40\;\textrm{dB}\right)\ldotp
% The stopband edge frequency \omega_s =0\ldotp 6\pi.
% The MATLAB function firhalfband is used for the FIR filter design.

clear all, close all

% Specifying the minimal stopband attenuation

as = 40;                              

% Computation of the peak ripple

delta = 10^(-as/20);                  

% Specifying the stopband edge frequency

fs = 0.60;                             

% Computation of the passband edge frequency

fp = 1-fs;                             

% Halfband filter design

h = firhalfband('minorder',fp,delta);  
figure 
stem(0:length(h)-1,h)
xlabel('n'), ylabel('h[n]')
title('Impulse response')
axis([0,22,-0.2,0.6])

% Halfband filter frequency response

[H,f] = freqz(h,1,512,2);
figure
plot(f,20*log10(abs(H))),grid
xlabel('\omega/\pi'), ylabel('Gain [dB]')
title('Gain response')
axis([0,1,-60,5])
disp('END OF EXAMPLE 7.2')

%
%% Example 7.3
% In this example , we design a halfband filter using the truncated impulse response window method. We take for the filter length N = 23, and the hamming window,

clear all,close all

% Specifying  the filter order

Nord = 22; 

% Specifying the window

win = hamming(23); 

% Halfband filter design

h = firhalfband(Nord,win); 
figure 
stem(0:length(h)-1,h)
xlabel('n'), ylabel('h[n]')
title('Impulse response')
axis([0,22,-0.2,0.6])

% Halfband filter frequency response

[H,f] = freqz(h,1,512,2);
figure
plot(f,20*log10(abs(H))),grid
xlabel('\omega/\pi'), ylabel('Gain [dB]')
title('Gain response')
axis([0,1,-60,5])
disp('END OF EXAMPLE 7.3')

%
%% Example 7.4
% This example demonstrates in MATLAB the process of decimation-by-2 based on the efficient structure shown bellow. The Figure shows polyphase memory saving implementation on the example of the linear-phase 10th order FIR halfband filter. The coefficient symmetry property is exploited.
%
% Efficient implementation of the linear-phase FIR halfbf and filter.
% (a) Polyphase configuration. (b) Implementation scheme for N=11.
%
% This program performs:
% Design of linear-phase FIR halfband filter of Nord = 22.
% Generates 257 samples of the original signal x
% Implements the decimation-by-2 according to the efficient scheme presented above. The implementation structure is extended since the filter ordeer in this example is 22 (instead of 11).The delay chain has 11 elements, and the number of multiplication constants is 6.
% Computes and plots the spectrum of the original signal and that of the decimated signal.
%

clear all, close all

% Designing 22nd order linear-phase FIR filter

Nord = 22; % Hafband filter order
h = firhalfband(22,0.4); % FIR filter design

% Generation of the input signal

F = [0,0.05,0.45,1]; A = [0,1,0.1,0.05];
x1 = fir2(256,F,A); x2 = sin(2*pi*(0:256)*0.4)/150;
x = x1 + x2;
[X,f] = freqz(x,1,512,2);  % Spectrum of the original signal

% Polyphase decomposition

e0 = h(1:2:(length(h)-1)/2);

% Setting the initial states

xk = zeros(size(1:(length(h)+1)/2));

% Decimation

rot = 0; % Initial swithch position
y = [];
for n=1:length(x)
xn = x(n);
   
if rot==0 ss=1;
    pro = xn*e0;
    xk = [pro+xk(1:length(xk)/2),fliplr(pro) + xk((length(xk)/2+1):length(xk))];
    yn = xk(length(xk));
    y = [yn,y];
    xk = [0,xk(1:length(xk)-1)];
else
end
if rot==1 ss=0;
    xk(length(xk)/2+1) = 0.5*xn + xk(length(xk)/2+1);
else
end
rot=ss;
end

% Displaying the results

figure 
subplot(2,1,1)
plot(f,abs(X(1:512))), xlabel('\omega/\pi'),ylabel('|X(e^{j\omega})|')
title('Spectrum of the input signal')
axis([0,1,0,1.1])
Y = freqz(y,1,512,2); % Spectrum of the decimated signal
subplot(2,1,2)
plot(f,abs(Y(1:512))), xlabel('\omega/\pi'),ylabel('|Y(e^{j\omega})|')
title('Spectrum of the decimated signal')
axis([0,1,0,0.6])
disp('END OF EXAMPLE 7.4')

%
%% 7.3 IIR Halfband Filter Design in MATLAB
%% Example 7.5
%
% Design an elliptic halfband filter that is specified with the stopband edge frequency of \omega_s =0\ldotp 57\pi and the minimal stopband attenuation a_s =45\;\textrm{dB}.
% For IIR halfband filter design we use the function halfbandiir which is included in this folder.
% The design procedure is performed in three steps:

clear all, close all

% Step 1: Adjusting the design parameters

as = 45
ap = 10*log10(1+1/(10^(as/10)-1))
fs = 0.57
fp = 1-fs

% Step 2: Computation of the minimal filtr order

[N,fp] = ellipord(fp,fs,ap,as)

% Step 3: Computation of the halfband filter parameters

[b,a,z,p,k] = halfbandiir(7,0.43);

% Halfband filter frequency response

[H,f]=freqz(b,a,1024,2);

% Setting the specifications

fspec=[0.57,0.57,1];
aspec=[2,-45,-45];

% Displaying the results

figure 
plot(f,20*log10(abs(H)))
hold on
plot(fspec,aspec,'r'), xlabel('\omega/\pi'), ylabel('Gain, dB')
axis([0,1,-80,2])
figure 
zplane(z,p)
disp('Notice that filter poles are located on the imaginary axis.')
disp('END OF EXAMPLE 7.5')

%
%% 7.4 IIR Halfband Filters with Approximately Linear Phase
%
% Remainder
% The main disadvantage of minimum phase IIR filters is their nonlinear phase response. A solution for the approximately linear-phase IIR halfband filter can be achieved when representing the transfer function H(z) as a sum of two subfilters, one of them a pure delay term,
% H\left(z\right)=\frac{1}{2}\left(A_0 \left(z\right)+z^{{-K}_1 } \right)
% where the subfilter A_0 \left(z\right) iz an all-pass function expressible as a product of the second-order and the fourth-order sections.
%
%% Example 7.6
%
% Design the approximately linear-phase IIR halfband filter specified by the following parameters:
% Given the passband/stopband edge frequencies: \omega_p =0\ldotp 4\pi, \omega_s =0\ldotp 6\pi.
% The all-pass branch A_0 \left(z\right) is of degree K_0 =10, and the degree of the delay branch is K_1 =9.
% Solution:
% The design method is based on the paper: H.W. Schussler, P. Stefen (2001), Recursive Halfband Filters, AEU Int. J. Electron. Commun., 55(6), 377-388.
% The paper was accompanied with the MATLAB program Hblipha, which were available at
% http://www-nt.e-technik.uni-erlangen.de/~hws/programs/halfbandfilters/ . The program Hblipha is included into this folder.
% Filter design

[b1,b2,a,aA,dP] = Hblipha(10,0.4,'Cheby');

% Frequency response

[H,f]=freqz(b1,a,1024,2);

% Group delay

[gd,f] = grpdelay(b1,aA,1024,2);

% Displaying the results

figure
plot(f,20*log10(abs(H)))
title('Gain response')
axis([0,1,-70,2]),xlabel('\omega/\pi'), ylabel('Gain, dB'),grid
figure
zplane(b1,aA)
title('pole-zero locations')
figure
plot(f,unwrap(angle(H))),xlabel('\omega/\pi'), ylabel('Phase, rad'),grid 
title('Phase response')
figure
plot(f,gd),,xlabel('\omega/\pi'), ylabel('Group delay, samples'),grid
title('Group delay')
disp('END OF EXAMPLE 7.6')
disp('                     END OF CHAPTER VII')

%

