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

%%            Output Info about this m-file

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

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

banner();

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

w1 = 0;       w2 = 0.20*pi; delta1 = 0.20; gain1 = 2.0;

w3 = 0.25*pi; w4 = 0.45*pi; delta2 = 0.05; gain2 = 0.0;

w5 = 0.55*pi; w6 = 0.70*pi; delta3 = 0.30; gain3 = 3.0;

w7 = 0.75*pi; w8 = pi;      delta4 = 0.10; gain4 = 1.0;

f = [w2, w3, w4, w5, w6, w7]/pi; m = [1, 0, 1, 1]; delta = [delta1, delta2, delta3, delta4];

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

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

fprintf('\n Results will be used by firpm function:\n');

N

%f0

%m0

%weights

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

weights = [delta3/delta1  delta3/delta2  1  delta3/delta4]

deltaH = max([delta1,delta2,delta3,delta4]); deltaL = min([delta1,delta2,delta3,delta4]);

%Dw1 = min((w3-w2), (w5-w4));

%Dw2 = min((w5-w4), (w7-w6));

%Dw  = min(Dw1, Dw2);

%M = ceil((-20*log10((delta1*delta2*delta3*delta4)^(1/4)) - 13) / (2.285*Dw) + 1)

f = [ 0   w2  w3  w4  w5  w6  w7  pi]/pi;

m = [ 2   2   0   0   3   3   1   1];

h = firpm(N+4, f, m, weights);               % even-number order

M = N+5

delta_w = 2*pi/1000;

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

w1i = floor(w1/delta_w)+1; w2i = floor(w2/delta_w)+1;

w3i = floor(w3/delta_w)+1; w4i = floor(w4/delta_w)+1;

w5i = floor(w5/delta_w)+1; w6i = floor(w6/delta_w)+1;

w7i = floor(w7/delta_w)+1; w8i = floor(w8/delta_w)+1;

Asd = -max(db(w3i:w4i))

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

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

l = 0:M-1;

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

%%                    Plot

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

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

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

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

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

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

set(gca,'YTickMode','manual','YTick',[-0.4:0.2:1.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,-34,-30,0]);

set(gca,'YTickLabelMode','manual','YTickLabel',['60';'34';'30';' 0']);

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

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

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

set(gca,'YTickMode','manual','YTick',[-0.05,0.05,0.9,1.1,1.8,2.2,2.7,3.3]);

delta_w = 2*pi/1000;

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

b1w = omega(w1i:w2i)/pi; b1e = (Hr(w1i:w2i)-m(1)); %b1e = (Hr(w1i:w2i)-m(1))*weights(1);

b2w = omega(w3i:w4i)/pi; b2e = (Hr(w3i:w4i)-m(3)); %b2e = (Hr(w3i:w4i)-m(3))*weights(2);

b3w = omega(w5i:w6i)/pi; b3e = (Hr(w5i:w6i)-m(5)); %b3e = (Hr(w5i:w6i)-m(5))*weights(3);

b4w = omega(w7i:w8i)/pi; b4e = (Hr(w7i:w8i)-m(7)); %b4e = (Hr(w7i:w8i)-m(7))*weights(4);

plot(b1w,b1e,b2w,b2e,b3w,b3e,b4w,b4e); grid on;

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

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

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

figure('NumberTitle', 'off', 'Name', 'Problem 7.35 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(Not to Scale)'); title('Amplitude Response');

set(gca,'YTickMode','manual','YTick',[-0.05,0.05,0.9,1.1,1.8,2.2,2.7,3.3]);

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

set(gca,'XTickMode','manual','XTick',[0,0.2,0.25,0.45,0.55,0.7,0.75,1]);