%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

%%            Output Info about this m-file

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

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

banner();

%% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

% bandpass

ws1 = 0.25*pi; wp1 = 0.4*pi; wp2 = 0.6*pi; ws2 = 0.75*pi;

As = 50; Rp = 1;

[delta1, delta2] = db2delta(Rp, As);

deltaH = max(delta1,delta2); deltaL = min(delta1,delta2);

f = [ws1, wp1, wp2, ws2]/pi; m = [0, 1, 0]; delta = [delta2, delta1, delta2];

[N, f, m, weights] = firpmord(f, m, delta);

N

h = firpm(N, f, m, weights);

[db, mag, pha, grd, w] = freqz_m(h, [1]);

delta_w = 2*pi/1000;

ws1i = floor(ws1/delta_w)+1; wp1i = floor(wp1/delta_w)+1;

ws2i = floor(ws2/delta_w)+1; wp2i = floor(wp2/delta_w)+1;

Asd = -max(db(1:1:ws1i))

N = N + 2            % Length M=25

h = firpm(N, f, m, weights);

[db, mag, pha, grd, w] = freqz_m(h, [1]);

Asd = -max(db(1:1:ws1i))

M = N + 1

l = 0:M-1;

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

%%                   Type-1 BPF

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

[Hr, ww, a, L] = Hr_Type1(h);

Rp = -(min(db(floor(wp1/delta_w)+1 :1: floor(wp2/delta_w))));     % Actual Passband Ripple

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

As = -round(max(db(floor(ws2/delta_w)+1 : 1 : 501)));        % Min Stopband attenuation

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

[delta1_db, delta2_db] = db2delta(Rp, As)

% Plot

figure('NumberTitle', 'off', 'Name', 'Problem 7.29 h(n), Parks-McClellan Method')

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

subplot(2,2,1); stem([0:M-1], h); axis([0 M-1 -0.5 0.5]); grid on;

xlabel('n'); ylabel('h(n)'); title('Actual Impulse Response, M=25');

subplot(2,2,2); plot(w/pi, db); axis([0 1 -90 10]); grid on;

set(gca,'YTickMode','manual','YTick',[-53,-9,0])

set(gca,'YTickLabelMode','manual','YTickLabel',['53';' 9';' 0']);

set(gca,'XTickMode','manual','XTick',[0,0.25,0.33,0.4,0.6,0.66,0.75,1]);

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

subplot(2,2,3); plot(ww/pi, Hr); axis([0, 1, -0.2, 1.2]); grid on;

xlabel('frequency in \pi nuits'); ylabel('Hr(w)'); title('Amplitude Response');

set(gca,'XTickMode','manual','XTick',[0,0.25,0.33,0.4,0.6,0.66,0.75,1])

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

subplot(2,2,4);

sb1w = ww(1:1:ws1i)/pi; sb1e = Hr(1:1:ws1i);

pbw = ww(wp1i:wp2i)/pi; pbe = Hr(wp1i:wp2i)-1;

sb2w = ww(ws2i:501)/pi; sb2e = Hr(ws2i:501);

plot(sb1w,sb1e,pbw,pbe*(delta2/delta1),sb2w,sb2e);   % weighted error

%plot(sb1w,sb1e,pbw,pbe,sb2w,sb2e);   %  error 

axis([0, 1, -deltaL, deltaL]); grid on;

xlabel('frequency in \pi units'); ylabel('Hr(w)');

title('Weighted Error');

%title('Error Response');

set(gca,'XTickMode','manual','XTick',f)

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

set(gca,'XGrid','on','YGrid','on')

figure('NumberTitle', 'off', 'Name', 'Problem 7.29a Parks-McClellan Method')

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

subplot(2,2,1); plot(w/pi, db); grid on; axis([0 2 -90 10]);

set(gca,'YTickMode','manual','YTick',[-53,-9,0])

set(gca,'YTickLabelMode','manual','YTickLabel',['53';' 9';' 0']);

set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,0.6,0.75,1,1.25,1.4,1.6,1.75,2]);

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

subplot(2,2,3); plot(w/pi, mag); grid on; %axis([0 1 -100 10]);

xlabel('frequency in \pi units'); ylabel('Absolute'); title('Magnitude Response in absolute');

set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,0.6,0.75,1,1.25,1.4,1.6,1.75,2]);

set(gca,'YTickMode','manual','YTick',[0,1.0]);

subplot(2,2,2); plot(w/pi, pha); grid on; %axis([0 1 -100 10]);

xlabel('frequency in \pi units'); ylabel('Rad'); title('Phase Response in Radians');

subplot(2,2,4); plot(w/pi, grd*pi/180);  grid on; %axis([0 1 -100 10]);

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

figure('NumberTitle', 'off', 'Name', 'Problem 7.29 AmpRes of h(n), Parks-McClellan Method')

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

plot(ww/pi, Hr); grid on; %axis([0 1 -100 10]);

xlabel('frequency in \pi units'); ylabel('Hr'); title('Amplitude Response');

set(gca,'YTickMode','manual','YTick',[-delta2_db ,0,delta2_db , 1-delta1_db, 1, 1+delta1_db]);

set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,0.6,0.75,1]);