%% ------------------------------------------------------------------------

%%            Output Info about this m-file

fprintf('\n***********************************************************\n');

fprintf('        <DSP using MATLAB> Problem 8.42 \n\n');

banner();

%% ------------------------------------------------------------------------

% Digital Filter Specifications:   Elliptic bandstop

wsbs = [0.40*pi 0.48*pi];             % digital stopband freq in rad

wpbs = [0.25*pi 0.75*pi];             % digital passband freq in rad

Rp = 1.0    % passband ripple in dB

As = 80      % stopband attenuation in dB                       

Ripple = 10 ^ (-Rp/20)           % passband ripple in absolute

Attn = 10 ^ (-As/20)             % stopband attenuation in absolute

% Calculation of Elliptic filter parameters:

[N, wn] = ellipord(wpbs/pi, wsbs/pi, Rp, As);

fprintf('\n  ********* Elliptic Digital Bandstop Filter Order is = %3.0f \n', 2*N)

% Digital Elliptic bandstop Filter Design:

[bbs, abs] = ellip(N, Rp, As, wn, 'stop');

[C, B, A] = dir2cas(bbs, abs)

% Calculation of Frequency Response:

[dbbs, magbs, phabs, grdbs, wwbs] = freqz_m(bbs, abs);

% ---------------------------------------------------------------

%    find Actual Passband Ripple and Min Stopband attenuation

% ---------------------------------------------------------------

delta_w = 2*pi/1000;

Rp_bs = -(min(dbbs(1:1:ceil(wpbs(1)/delta_w+1))));                           % Actual Passband Ripple

fprintf('\nActual Passband Ripple is %.4f dB.\n', Rp_bs);

As_bs = -round(max(dbbs(ceil(wsbs(1)/delta_w)+1:1:ceil(wsbs(2)/delta_w)+1)));   % Min Stopband attenuation

fprintf('\nMin Stopband attenuation is %.4f dB.\n\n', As_bs);

%% -----------------------------------------------------------------

%%                             Plot

%% -----------------------------------------------------------------  

figure('NumberTitle', 'off', 'Name', 'Problem 8.42 Elliptic bs by ellip function')

set(gcf,'Color','white');

M = 1;                          % Omega max

subplot(2,2,1); plot(wwbs/pi, magbs); axis([0, M, 0, 1.2]); grid on;

xlabel('Digital frequency in \pi units'); ylabel('|H|'); title('Magnitude Response');

set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.25, 0.40, 0.48, 0.75, M]);

set(gca, 'YTickMode', 'manual', 'YTick', [0, 0.01, 0.8913, 1]);

subplot(2,2,2); plot(wwbs/pi, dbbs); axis([0, M, -120, 2]); grid on;

xlabel('Digital frequency in \pi units'); ylabel('Decibels'); title('Magnitude in dB');

set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.25, 0.40, 0.48, 0.75, M]);

set(gca, 'YTickMode', 'manual', 'YTick', [ -80, -40, 0]);

set(gca,'YTickLabelMode','manual','YTickLabel',['80'; '40';' 0']);

subplot(2,2,3); plot(wwbs/pi, phabs/pi); axis([0, M, -1.1, 1.1]); grid on;

xlabel('Digital frequency in \pi nuits'); ylabel('radians in \pi units'); title('Phase Response');

set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.25, 0.40, 0.48, 0.75, M]);

set(gca, 'YTickMode', 'manual', 'YTick', [-1:0.5:1]);

subplot(2,2,4); plot(wwbs/pi, grdbs); axis([0, M, 0, 50]); grid on;

xlabel('Digital frequency in \pi units'); ylabel('Samples'); title('Group Delay');

set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.25, 0.40, 0.48, 0.75, M]);

set(gca, 'YTickMode', 'manual', 'YTick', [0:20:50]);

% ------------------------------------------------------------

%                 PART 2

% ------------------------------------------------------------

% Discrete time signal

Ts = 1;                              % sample intevals

n1_start = 0; n1_end = 200;

      n1 = [n1_start:n1_end];        % [0:200]

xn1 = sin(0.44*pi*n1);               % digital signal

% ----------------------------

%       DTFT of xn1

% ----------------------------

M = 500;

[X1, w] = dtft1(xn1, n1, M);

%magX1 = abs(X1);

angX1  = angle(X1);  realX1  = real(X1);  imagX1  = imag(X1);

magX1 = sqrt(realX1.^2 + imagX1.^2);

%% --------------------------------------------------------------------

%%              START X(w)'s  mag ang real imag

%% --------------------------------------------------------------------

figure('NumberTitle', 'off', 'Name', 'Problem 8.42 X1 DTFT');

set(gcf,'Color','white');

subplot(2,1,1); plot(w/pi,magX1); grid on;  %axis([-1,1,0,1.05]);

title('Magnitude Response');

xlabel('digital frequency in \pi units'); ylabel('Magnitude  |H|');

set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.2, 0.44, 0.6, 0.8, 1.0, 1.2, 1.4, 1.56, 1.8, 2]);

subplot(2,1,2); plot(w/pi, angX1/pi); grid on;  %axis([-1,1,-1.05,1.05]);

title('Phase Response');

xlabel('digital frequency in \pi units'); ylabel('Radians/\pi');

figure('NumberTitle', 'off', 'Name', 'Problem 8.42 X1 DTFT');

set(gcf,'Color','white');

subplot(2,1,1); plot(w/pi, realX1); grid on;

title('Real Part');

xlabel('digital frequency in \pi units'); ylabel('Real');

set(gca, 'XTickMode', 'manual', 'XTick', [0, 0.2, 0.44, 0.6, 0.8, 1.0, 1.2, 1.4, 1.56, 1.8, 2]);

subplot(2,1,2); plot(w/pi, imagX1); grid on;

title('Imaginary Part');

xlabel('digital frequency in \pi units'); ylabel('Imaginary');

%% -------------------------------------------------------------------

%%             END X's  mag ang real imag

%% -------------------------------------------------------------------

% ------------------------------------------------------------

%                 PART 3

% ------------------------------------------------------------

yn1 = filter(bbs, abs, xn1);

n2 = n1;

% ----------------------------

%       DTFT of yn1

% ----------------------------

M = 500;

[Y1, w] = dtft1(yn1, n2, M);

%magY1  = abs(Y1);

angY1  = angle(Y1);  realY1  = real(Y1);  imagY1  = imag(Y1);

magY1 = sqrt(realY1.^2 + imagY1.^2);

%% --------------------------------------------------------------------

%%              START Y1(w)'s  mag ang real imag

%% --------------------------------------------------------------------

figure('NumberTitle', 'off', 'Name', 'Problem 8.42 Y1 DTFT');

set(gcf,'Color','white');

subplot(2,1,1); plot(w/pi,magY1); grid on;  %axis([-1,1,0,1.05]);

title('Magnitude Response');

xlabel('digital frequency in \pi units'); ylabel('Magnitude  |H|');

subplot(2,1,2); plot(w/pi, angY1/pi); grid on;  %axis([-1,1,-1.05,1.05]);

title('Phase Response');

xlabel('digital frequency in \pi units'); ylabel('Radians/\pi');

figure('NumberTitle', 'off', 'Name', 'Problem 8.42 Y1 DTFT');

set(gcf,'Color','white');

subplot(2,1,1); plot(w/pi, realY1); grid on;

title('Real Part');

xlabel('digital frequency in \pi units'); ylabel('Real');

subplot(2,1,2); plot(w/pi, imagY1); grid on;

title('Imaginary Part');

xlabel('digital frequency in \pi units'); ylabel('Imaginary');

%% -------------------------------------------------------------------

%%             END Y1's  mag ang real imag

%% -------------------------------------------------------------------

figure('NumberTitle', 'off', 'Name', 'Problem 8.42 x(n) and y(n)')

set(gcf,'Color','white');

subplot(2,1,1); stem(n1, xn1);

xlabel('n'); ylabel('x(n)');

title('xn sequence');  grid on;

subplot(2,1,2); stem(n1, yn1);

xlabel('n'); ylabel('y(n)');

title('yn sequence');  grid on;