
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % .m script of web demo "IIR Polyphase Filter"
    % Written for GNU Octave
    %
    % INSTITUTE OF TELECOMMUNICATIONS
    % University of Stuttgart
    % www.inue-uni-stuttgart.de
    % author: Sagar Hukkire
    % date: 14.01.2016
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % This example script just prints the values of the given parameters.
    % Make sure that the function prototype matches the prototype defined
    % in the XML config file. (Parameter names are arbitrary but should be the same
    % for better readability)

    function [] =IIR_Polyphase(samplingFrequency,DigitalCutoff,order,display_option,filename)



    % define formats
    landscape = "-S930,350";
    portrait  = "-S640,480";

    % select format
    output_format = landscape;

    % Create an invisible figure.
    fh = figure(1);
    set(fh, "visible", "off");

    % define style
    set(0, "defaultlinelinewidth", 2);
    set(0, "defaultaxesfontsize", 8);
    set(0, "defaultaxesfontname", "Helvetica");
    set(0, "defaulttextfontsize", 8);
    set(0, "defaulttextfontname", "Helvetica");






    %Analog_to_digital_filter conversion starts


    t=1:order;
    analogPoles=exp((1j*(2*t+order-1)*pi)./(2*order));
    den=abs(poly(analogPoles));
    analogZeros=ones(1,order)*1e10;
    fs=samplingFrequency;
    fp=DigitalCutoff;
    u0=1;
    k=1;


    k=tan((2*pi*fp)/(2*fs));
    y=u0/k;

    for i=1:order
        digitalPoles(:,i)=(1+(analogPoles(i)/y))./(1-(analogPoles(i)/y));
        digitalPolesCoeff=poly(digitalPoles(:));
        digitalPolesCoeff=abs(digitalPolesCoeff).*sign(real(digitalPolesCoeff));%digital_poles_coefficient
    end

    for i=1:order
        digitalZeros(:,i)=(1+(analogZeros(i)/y))./(1-(analogZeros(i)/y));
        digitalZerosCoeff=abs(poly(digitalZeros(:)));
    end
    [tfdigitalDesign,w]=freqz(digitalZerosCoeff,digitalPolesCoeff);%digital_zeros

    tfdigitalZero=freqz(digitalZerosCoeff,1);
    tfdigitalpole=freqz(digitalPolesCoeff,1);
    if(display_option=='Dt')


        subplot(2,1,1);

        [hax]=plotyy(w/pi,abs(tfdigitalZero),w/pi,angle(tfdigitalZero));
        

        xlabel('\omega / \pi');

        ylabel(hax(1),'Magnitude');

        ylabel(hax(2),'Phase');

        legend('Numerator polynomial frequency response');

        subplot(2,1,2);


        [hax]=plotyy(w/pi,abs(tfdigitalpole),w/pi,angle(tfdigitalpole));
        


        xlabel('\omega / \pi');

        ylabel(hax(1),'Magnitude');

        ylabel(hax(2),'Phase');

        legend('Denominator polynomial frequency response');




    end
    %Analog_to_digital_filter conversion ends



    %IIR Polyphase code starts

    %Zero padding for numerator starts

    lenZeros=length(digitalZerosCoeff);
    M=2;
    Qutiont=fix(lenZeros/M);

    temp=lenZeros-(Qutiont.*M);

    padvalue=M-temp;

    numZeroPadd=[digitalZerosCoeff zeros(1,padvalue)];

    %Zero padding for numerator ends

    %Zero padding for denominator starts

    lenPoles=length(digitalPolesCoeff);
    M=2;
    Qutiont1=fix(lenPoles/M);

    temp1=lenZeros-(Qutiont1.*M);

    padvalue1=M-temp1;

    denPolePadd=[-digitalPolesCoeff zeros(1,padvalue1)];

    %Zero padding for denominator ends




    E=reshape(numZeroPadd,M,length(numZeroPadd)/M);%subband filters of numerator
    H=reshape(denPolePadd,M,length(denPolePadd)/M);%subband filters of denominator



    [new1,w]=freqz(E(1,:),1,1024,'whole');%frequency response of suband filter

    [new2,w]=freqz(E(2,:),1,1024,'whole');

    [new3,w]=freqz(H(1,:),1,1024,'whole');

    [new4,w]=freqz(H(2,:),1,1024,'whole');

    inv_z=exp(-(j*w/2));

    if(display_option=='Sd')

        subplot(2,2,1);

        [hax]=plotyy((0.5*w)/pi,abs(new1),(0.5*w)/pi,angle(new1));
        

        xlabel('\omega / \pi');

        ylabel(hax(1),'Magnitude');

        legend('Numerator Polynomial 1');

        



        subplot(2,2,2);
        [hax]=plotyy((0.5*w)/pi,abs(new2.*inv_z),(0.5*w)/pi,angle(new2.*inv_z));
        

        xlabel('\omega / \pi');

        

        ylabel(hax(2),'Phase');
        legend('Numerator Polynomial 2');

        subplot(2,2,3);
        [hax]=plotyy((0.5*w)/pi,abs(new3),(0.5*w)/pi,angle(new3));
        
        xlabel('\omega / \pi');
        ylabel(hax(1),'Magnitude');
        legend('denominator Polynomial 1');
        


        subplot(2,2,4);
        [hax]=plotyy((0.5*w)/pi,abs(new4.*inv_z),(0.5*w)/pi,angle(new4.*inv_z));
        
        xlabel('\omega / \pi');
        
        ylabel(hax(2),'Phase');
        legend('denominator Polynomial 2');
    end



    new5=new1+(new2.*inv_z);
    new6=new3+(new4.*inv_z);
    new7=new5./new6;

    if(display_option=='Pl')
        plot((0.5*w)/pi,abs(new7)) %polyphase response
        hold on

        plot([0, 2*fp/fs], [0.707*max(abs(new7)), 0.707*max(abs(new7))], 'r');
        plot([2*fp/fs, 2*fp/fs], [0, 0.707*max(abs(new7))], 'r');

        xlabel('\omega / \pi');
        ylabel('|H(\omega)|');
        grid on
        grid minor
    end








    % Create an image file. This png image is displayed on the webpage. Its size
    % must be 930x350.
    print(filename, "-dpng", output_format);

    end

