% Parks MacClellan optimization alggorithm
pkg load signal
order=80;
f=[0 0.465 0.535 1.0]; % frequências qie limitam as bandas
m=[1 1 0 1]; % Ganhos do filtro às frequências que limitam as bandas
w=[1 1];% Peso de cada banda
h=remez(order,f,m,w); % em MTALAB h=firpm(order,f,m,w)
close all

% H(frequencia normalizada)
figure(1)
frange=pi*[-1:1/256:1-1/256];
[H,W]=freqz(h,1,frange);
plot(W/pi,abs(H))
xlabel('\omega/\pi')
ylabel('GAIN')
grid

% % H(frequencia em Khz)
figure(2)
fs=8;
f=W/(2*pi)*fs;
plot(f(257:end),abs(H(257:end)))
xlabel('f(kHz)')
ylabel('GAIN')
grid

% H(frequência normalizada de 0 a pi)
figure(3)
plot(W(257:end)/pi,abs(H(257:end)))
xlabel('\omega/\pi')
ylabel('GAIN')
grid


% H(z^2) frequência normalizada
h2=zeros(2*order+1,1);

for k=1:order+1
  h2(2*k-1)=h(k);
end

[H2 W2]=freqz(h2,1,frange);
figure(4)
plot(W2/pi,abs(H2))
xlabel('\omega/\pi')
ylabel('|H(z^2)|')
grid

% H(-z) frequência normalizada
h3=zeros(size(h));
for k=1:order+1
  h3(k)=h(k)*(-1)^(k-1)
end
figure(5)
[H3 W3]=freqz(h3,1,frange);
plot(W3/pi,abs(H3))
xlabel('\omega/\pi')
ylabel('|H(-z)|')
grid