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

%%            Output Info about this m-file

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

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

banner();

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

% bandpass, PM method

ws1 = 0.2*pi; wp1 = 0.35*pi; wp2 = 0.55*pi; ws2 = 0.75*pi;

Rp = 0.25; As = 40;

[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, f0, m0, weights] = firpmord(f, m, delta);

fprintf('\n-------------------------------------------\n');

fprintf('\n Results by firpmord function:\n');

N

%f0

%m0

%weights

fprintf('\n-------------------------------------------\n');

weights = [delta1/delta2   1   delta1/delta2]

f = [ 0     ws1    wp1    wp2   ws2   pi]/pi;

m = [ 0     0      1      1     0     0];

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

M = N + 1

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

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

M = N + 3

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

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

[Hr,omega,P,L] = ampl_res(h);

l = 0:M-1;

delta_w = 2*pi/1000; sp_edge1 = ws1/delta_w+1; sp_edge2 = ws2/delta_w+1;

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

%%                    Plot

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

figure('NumberTitle', 'off', 'Name', 'Problem 7.38')

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

subplot(2,2,1); stem(l, h); axis([-1, M, -0.4, 0.45]); grid on;

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

set(gca,'XTickMode','manual','XTick',[0,(M-1)/2,M-1]);

set(gca,'YTickMode','manual','YTick',[-0.4:0.1:0.4]);

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

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

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

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

set(gca,'YTickLabelMode','manual','YTickLabels',['60';'40';' 0']);

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

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

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

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

subplot(2,2,4); axis([0, 1, -deltaH, deltaH]);

sb1w = omega(1:1:ws1i)/pi;  sb1e = Hr(1:1:ws1i)-m(1);  %sb1e = (Hr(1:1:ws1i)-m(1))*weights(1);

pbw  = omega(wp1i:wp2i)/pi; pbe = Hr(wp1i:wp2i)-m(3);  %pbe  = (Hr(wp1i:wp2i)-m(3))*weights(2);

sb2w = omega(ws2i:501)/pi;  sb2e = Hr(ws2i:501)-m(5);  %sb2e = (Hr(ws2i:501)-m(5))*weights(3);   

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

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

xlabel('frequency in \pi units'); ylabel('Hr(\omega)'); title('Error Response'); %title('Weighted Error');  

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

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

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

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

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

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

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

set(gca,'YTickMode','manual','YTick',[-delta2,delta2, 1 - delta1, 1 + delta1]);

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