% racunanje koeficijenata analizirajuceg i rekonstrukcijskog filtarskog
% sloga za pseudokvadraturni filtarski slog s gotovo potpunom
% rekonstrukcijom (eng. pseudo QMF)
% ulazni parametar M predstavlja broj pojaseva filtarskog sloga
% ulazni parametar N predstavlja red pojedinog filtra u filtarskom slogu
% NP prototip se konstruira uporabom Parks-McClellan algoritma za
% konstrukciju FIR filtra

function [h,f]=pseudoqmf(M,N)

% definiranje ulaznih parametara NP prototipa
a=[1 1 0 0]; % definiranje zeljenih amplituda
weight=[10 1]; % definiranje tezinskih faktora
fakt=1; % inicijalni faktor za odredjivane granice pojasa propustanja
eta=.005; % zeljena max. valovitost u podrucju propustanja
step=1; % pocetni korak optimizacije
dir=1; % pocetni smjer optimizacije
err=1; % inicijalna pogreska 
prev_err=50; % inicijalna prethodna pogreska
brojac=0; % brojac u optimizacijskoj petlji

while ((err>eta))
    if (brojac>100) % max. broj iteracija ogranicen na 100
        break
    end
    brojac=brojac+1;
    % provjera je li iteracija povecala pogresku i korekcija parametara ako
    % jest povecala pogresku
    if (err>=prev_err)
        step=step/2;
        dir=-dir;
    end

    % promjena faktora kojim se definira granica pojasa propustanja NP
    % prototipa
    fakt=fakt+step*dir;
    
    f=[0 1/(fakt*M) 1/M 1];
        
    % kreiranje prototipa
    b=firpm(N,f,a,weight);
    
    % konstrukcija sloga za analizu
    for k=0:M-1
        for j=0:N
            h(k+1,j+1)=2*b(j+1)*cos((2*k+1)*pi/(2*M)*(j-(N-1)/2)+(-1)^k*(pi/4));
        end
        h_d(k+1,:)=downsample(h(k+1,:),2);
        [hf(k+1,:),w(k+1,:)]=freqz(h(k+1,:));
    end
    
    % konstrukcija sloga za sintezu
    for k=0:M-1
        for j=0:N
            f(k+1,j+1)=2*b(j+1)*cos((2*k+1)*pi/(2*M)*(j-(N-1)/2)-(-1)^k*(pi/4));
        end
        [ff(k+1,:),wf(k+1,:)]=freqz(f(k+1,:));
    end
    
    %kaskadiranje sloga za analizu i sloga za sintezu
    izlaz_s=zeros(1,2*(N+1)-1);
    for k=1:M
        izlaz(k,:)=conv(h(k,:),f(k,:));
        izlaz_s=izlaz_s+(izlaz(k,:));
    end
    izlaz_h=20*log10(abs(fft(izlaz_s)));
    % odredjivanje pogreske filtra
    min_err=min(abs(izlaz_h));
    max_err=max(abs(izlaz_h));
    prev_err=err;
    err=max(min_err,max_err);
end