% cosine_desing - Funkcija odreduje 'dizajn' prototipnog filtra za kosinusni
%                 dizajn filtarskih slogova koristeci vremenske otvore.
%               - koristi se closed-form window metoda
%
% Poziva se:
%   [hp,h,N] = codesing(M,Aa,ro,izbor1,izbor2);
%
% gdje je:
%   M      - broj kanala (pojaseva; filtarska banka)
%   Aa     - atenuacija u pojasu gušenja (stopband attenuation)-unijeti
%            sa pozitivnim predznakom
%   ro     - parametar dizajna; roll-off faktor-pokazuje preklapanje izmedu 
%            pojaseva filtara
%   pobuda - pobudni signal
%
%
%  Primjer poziva:
%   [hp,h,N] = codesign(10,40,1.05,2,2);

function [hp,xr,N] = codesign(M,Aa,ro,x,izbor1,izbor2)

% -------------------------------------------------------------------------
%                   dizajn CMFBs-a koristeci window metodu
% -------------------------------------------------------------------------

x_save = x;
% matrica koja sadrzi koeficijente a1,b1,c1,a2,b2,c2 za ralicite odabire
% parametra Aa
mat = [ 50  3.096e-2 -6.345e-1 1.015     8.307e-2 -8.016e-1 1.130;
        60  2.150e-3 -6.286e-1 1.002    -1.543e-2 -5.931e-1 9.838e-1;
        70 -1.173e-3 -6.502e-1 9.998e-1 -1.909e-2 -6.043e-1 9.711e-1;
        80  3.107e-4 -6.735e-1 1.000    -3.959e-4 -6.738e-1 1.002;
        90 -1.839e-3 -6.880e-1 9.994e-1 -9.398e-3 -6.680e-1 9.866e-1;
       100  1.504e-3 -7.074e-1 1.001     7.341e-3 -7.237e-1 1.011;
       110 -7.235e-4 -7.167e-1 9.996e-1  1.756e-2 -7.601e-1 1.025;
       120  2.536e-3 -7.329e-1 1.001    -9.054e-3 -7.085e-1 9.882e-1;
       130  4.979e-4 -7.398e-1 1.000     6.255e-3 -7.547e-1 1.010;
       140  2.707e-4 -7.484e-1 1.000    -4.075e-3 -7.403e-1 9.964e-1;
       150 -4.951e-4 -7.553e-1 9.999e-1 -1.072e-2 -7.322e-1 9.870e-1];
   
Aa_new = 10*(round(Aa/10));

% za dani Aa racunaju se potrebni koeficijetni koji su rezultat
% optimizacije
for i = 1:length(mat)
    if (mat(i,1) == Aa_new)
        a1 = mat(i,2);
        b1 = mat(i,3);
        c1 = mat(i,4);
        a2 = mat(i,5);
        b2 = mat(i,6);
        c2 = mat(i,7);
        break;
    end
end

% racunanje granicne frekvencije podrucja gusenja
wa = ((1+ro)*pi)/(2*M);  

% racunanje parametra D (empirijska formula; rezultat optimizacije)
if ((ro < 1.0) & (ro >= 0.5))
	D = a1*(ro^2) + b1*ro + c1;
end
if ((ro <= 1.5) & (ro >= 1.0))
	D = a2*(ro^2) + b2*ro+c2;
end

% racunanje granicne frekvencije podrucja propustanja (optimalni wp)
wp = D*pi/(2*M);


% racunanje parametra beta
if Aa > 50,
   beta = 0.1102*(Aa-8.7);
elseif((Aa <= 50) & (Aa >= 21))
    beta = 0.5842*(Aa-21)^0.4+0.07886*(Aa-21);
else 
    beta = 0;
end

delta_w = wa-wp;

% racunanje reda prototipnog filtra (NP filtar) koji se kosinusno modulira
N = floor((Aa-7.95)/(14.36*(delta_w/(2*pi))));

% racunanje centralne frekvencije
wc = (wp + wa)/2;

% racunanje prototipnog filtra-Kaiserov otvor
hp = fir1(N,wc/(pi),kaiser(N+1,beta));

% racunanje analizirajucih i rekonstrukcijiskih filtara
h = zeros(M,N+1);
f = zeros(M,N+1);
for k = 0:M-1
    for n = 0:N
        h(k+1,n+1) = 2*hp(n+1)*cos(pi/M*(k+0.5)*(n-N/2)+(-1)^k*pi/4);
        f(k+1,N+1-n) = h(k+1,n+1);    
    end
 end


% racunanje dimenzija matrice za pohranu podataka 
h_temp = conv(h(5,:),x);
tn = length(h_temp);
hint = zeros(M,tn);


if(izbor2 == 1)
    % racunanje konvolucija, decimacija i interpolacija
    for i = 1:M
        h_temp = conv(h(i,:),x);
        for s = 1:length(h_temp)
            if((mod(s,M)==0)||(s==1))
                h_tempint(s) = h_temp(s);
            else
                h_tempint(s) = 0;
            end
            hint(i,s) = h_tempint(s); 
        end
    end
else
    % racunanje konvolucija bez decimacije i interpolacije
    for i = 1:M
        h_temp = conv(h(i,:),x);
        for s = 1:length(h_temp)
            hint(i,s) = h_temp(s); 
        end
    end
end 
   
% rekonstruiranni signal-konvolucija
xr = conv(hint(1,:),f(1,:));
for i = 2:M
xr = xr+conv(hint(i,:),f(i,:));    
end

% racun funkcije dzistorzije
To=zeros(M,2*N+1);
for m=0:M-1
     for k=0:M-1
        To(m+1,:)=To(m+1,:)+conv(h(k+1,:),f(k+1,:));
     end
end


% prikaz funkcije distorzije
figure(10);
N1 = 512;
[H,w]=freqz(To(1,:),1,N1*M/2);
plot(w(1:N1/2)/pi,abs(H(1:N1/2)),'m'), grid;