function f = stepresp(x, a, b, c,d,M1)%STEPRESP computes the step response of an analog filter with any number of pole pairs.%%	This function takes four arguments:%	  x is time%	  a are the real values of the poles%	  b are imaginary values of the poles%	  c are the real values of the zeros%	  c are imaginary values of the zeros%	  M1 scales the step responseexpr1 = b(1).^(-1);expr2 = a(2).^2 + b(2).^2;expr3 = a(3).^2 + b(3).^2;expr4 = a(4).^2 + b(4).^2;expr5 = (a(1) - a(2)).^2 + (b(1) - b(2)).^2;expr6 = (a(1) - a(2)).^2 + (b(1) + b(2)).^2;expr7 = (a(1) - a(3)).^2 + (b(1) - b(3)).^2;expr8 = (a(1) - a(3)).^2 + (b(1) + b(3)).^2;expr9 = (a(1) - a(4)).^2 + (b(1) - b(4)).^2;expr10 = (a(1) - a(4)).^2 + (b(1) + b(4)).^2;expr11 = -1./2;expr12 = (c(1).^2 + d(1).^2).^(-1);expr13 = (c(2).^2 + d(2).^2).^(-1);expr14 = (c(3).^2 + d(3).^2).^(-1);expr15 = (c(4).^2 + d(4).^2).^(-1);expr16 = (a(1) - c(1)).^2 + (b(1) - d(1)).^2;expr17 = (a(1) - c(1)).^2 + (b(1) + d(1)).^2;expr18 = (a(1) - c(2)).^2 + (b(1) - d(2)).^2;expr19 = (a(1) - c(2)).^2 + (b(1) + d(2)).^2;expr20 = (a(1) - c(3)).^2 + (b(1) - d(3)).^2;expr21 = (a(1) - c(3)).^2 + (b(1) + d(3)).^2;expr22 = (a(1) - c(4)).^2 + (b(1) - d(4)).^2;expr23 = (a(1) - c(4)).^2 + (b(1) + d(4)).^2;expr24 = 1./2;expr25 = b(1);expr26 = angle(1 + (i.*(b(1) - b(2)))./(a(1) - a(2)));expr27 = angle(1 + (i.*(b(1) + b(2)))./(a(1) - a(2)));expr28 = angle(1 + (i.*(b(1) - b(3)))./(a(1) - a(3)));expr29 = angle(1 + (i.*(b(1) + b(3)))./(a(1) - a(3)));expr30 = angle(1 + (i.*(b(1) - b(4)))./(a(1) - a(4)));expr31 = angle(1 + (i.*(b(1) + b(4)))./(a(1) - a(4)));expr32 = -angle(1 + (i.*(b(1) - d(1)))./(a(1) - c(1)));expr33 = -angle(1 + (i.*(b(1) + d(1)))./(a(1) - c(1)));expr34 = -angle(1 + (i.*(b(1) - d(2)))./(a(1) - c(2)));expr35 = -angle(1 + (i.*(b(1) + d(2)))./(a(1) - c(2)));expr36 = -angle(1 + (i.*(b(1) - d(3)))./(a(1) - c(3)));expr37 = -angle(1 + (i.*(b(1) + d(3)))./(a(1) - c(3)));expr38 = -angle(1 + (i.*(b(1) - d(4)))./(a(1) - c(4)));expr39 = -angle(1 + (i.*(b(1) + d(4)))./(a(1) - c(4)));expr40 = a(1);expr41 = exp(x.*a(1));expr42 = cos(x.*b(1));expr43 = a(1).^2 + b(1).^2;expr44 = sin(x.*b(1));expr45 = b(2).^(-1);expr46 = (-a(1) + a(2)).^2 + (-b(1) + b(2)).^2;expr47 = (-a(1) + a(2)).^2 + (b(1) + b(2)).^2;expr48 = (a(2) - a(3)).^2 + (b(2) - b(3)).^2;expr49 = (a(2) - a(3)).^2 + (b(2) + b(3)).^2;expr50 = (a(2) - a(4)).^2 + (b(2) - b(4)).^2;expr51 = (a(2) - a(4)).^2 + (b(2) + b(4)).^2;expr52 = (a(2) - c(1)).^2 + (b(2) - d(1)).^2;expr53 = (a(2) - c(1)).^2 + (b(2) + d(1)).^2;expr54 = (a(2) - c(2)).^2 + (b(2) - d(2)).^2;expr55 = (a(2) - c(2)).^2 + (b(2) + d(2)).^2;expr56 = (a(2) - c(3)).^2 + (b(2) - d(3)).^2;expr57 = (a(2) - c(3)).^2 + (b(2) + d(3)).^2;expr58 = (a(2) - c(4)).^2 + (b(2) - d(4)).^2;expr59 = (a(2) - c(4)).^2 + (b(2) + d(4)).^2;expr60 = b(2);expr61 = angle(1 + (i.*(-b(1) + b(2)))./(-a(1) + a(2)));expr62 = angle(1 + (i.*(b(1) + b(2)))./(-a(1) + a(2)));expr63 = angle(1 + (i.*(b(2) - b(3)))./(a(2) - a(3)));expr64 = angle(1 + (i.*(b(2) + b(3)))./(a(2) - a(3)));expr65 = angle(1 + (i.*(b(2) - b(4)))./(a(2) - a(4)));expr66 = angle(1 + (i.*(b(2) + b(4)))./(a(2) - a(4)));expr67 = -angle(1 + (i.*(b(2) - d(1)))./(a(2) - c(1)));expr68 = -angle(1 + (i.*(b(2) + d(1)))./(a(2) - c(1)));expr69 = -angle(1 + (i.*(b(2) - d(2)))./(a(2) - c(2)));expr70 = -angle(1 + (i.*(b(2) + d(2)))./(a(2) - c(2)));expr71 = -angle(1 + (i.*(b(2) - d(3)))./(a(2) - c(3)));expr72 = -angle(1 + (i.*(b(2) + d(3)))./(a(2) - c(3)));expr73 = -angle(1 + (i.*(b(2) - d(4)))./(a(2) - c(4)));expr74 = -angle(1 + (i.*(b(2) + d(4)))./(a(2) - c(4)));expr75 = a(2);expr76 = exp(x.*a(2));expr77 = cos(x.*b(2));expr78 = sin(x.*b(2));expr79 = b(3).^(-1);expr80 = (-a(1) + a(3)).^2 + (-b(1) + b(3)).^2;expr81 = (-a(1) + a(3)).^2 + (b(1) + b(3)).^2;expr82 = (-a(2) + a(3)).^2 + (-b(2) + b(3)).^2;expr83 = (-a(2) + a(3)).^2 + (b(2) + b(3)).^2;expr84 = (a(3) - a(4)).^2 + (b(3) - b(4)).^2;expr85 = (a(3) - a(4)).^2 + (b(3) + b(4)).^2;expr86 = (a(3) - c(1)).^2 + (b(3) - d(1)).^2;expr87 = (a(3) - c(1)).^2 + (b(3) + d(1)).^2;expr88 = (a(3) - c(2)).^2 + (b(3) - d(2)).^2;expr89 = (a(3) - c(2)).^2 + (b(3) + d(2)).^2;expr90 = (a(3) - c(3)).^2 + (b(3) - d(3)).^2;expr91 = (a(3) - c(3)).^2 + (b(3) + d(3)).^2;expr92 = (a(3) - c(4)).^2 + (b(3) - d(4)).^2;expr93 = (a(3) - c(4)).^2 + (b(3) + d(4)).^2;expr94 = b(3);expr95 = angle(1 + (i.*(-b(1) + b(3)))./(-a(1) + a(3)));expr96 = angle(1 + (i.*(b(1) + b(3)))./(-a(1) + a(3)));expr97 = angle(1 + (i.*(-b(2) + b(3)))./(-a(2) + a(3)));expr98 = angle(1 + (i.*(b(2) + b(3)))./(-a(2) + a(3)));expr99 = angle(1 + (i.*(b(3) - b(4)))./(a(3) - a(4)));expr100 = angle(1 + (i.*(b(3) + b(4)))./(a(3) - a(4)));expr101 = -angle(1 + (i.*(b(3) - d(1)))./(a(3) - c(1)));expr102 = -angle(1 + (i.*(b(3) + d(1)))./(a(3) - c(1)));expr103 = -angle(1 + (i.*(b(3) - d(2)))./(a(3) - c(2)));expr104 = -angle(1 + (i.*(b(3) + d(2)))./(a(3) - c(2)));expr105 = -angle(1 + (i.*(b(3) - d(3)))./(a(3) - c(3)));expr106 = -angle(1 + (i.*(b(3) + d(3)))./(a(3) - c(3)));expr107 = -angle(1 + (i.*(b(3) - d(4)))./(a(3) - c(4)));expr108 = -angle(1 + (i.*(b(3) + d(4)))./(a(3) - c(4)));expr109 = a(3);expr110 = exp(x.*a(3));expr111 = cos(x.*b(3));expr112 = sin(x.*b(3));expr113 = b(4).^(-1);expr114 = (-a(1) + a(4)).^2 + (-b(1) + b(4)).^2;expr115 = (-a(1) + a(4)).^2 + (b(1) + b(4)).^2;expr116 = (-a(2) + a(4)).^2 + (-b(2) + b(4)).^2;expr117 = (-a(2) + a(4)).^2 + (b(2) + b(4)).^2;expr118 = (-a(3) + a(4)).^2 + (-b(3) + b(4)).^2;expr119 = (-a(3) + a(4)).^2 + (b(3) + b(4)).^2;expr120 = (a(4) - c(1)).^2 + (b(4) - d(1)).^2;expr121 = (a(4) - c(1)).^2 + (b(4) + d(1)).^2;expr122 = (a(4) - c(2)).^2 + (b(4) - d(2)).^2;expr123 = (a(4) - c(2)).^2 + (b(4) + d(2)).^2;expr124 = (a(4) - c(3)).^2 + (b(4) - d(3)).^2;expr125 = (a(4) - c(3)).^2 + (b(4) + d(3)).^2;expr126 = (a(4) - c(4)).^2 + (b(4) - d(4)).^2;expr127 = (a(4) - c(4)).^2 + (b(4) + d(4)).^2;expr128 = b(4);expr129 = angle(1 + (i.*(-b(1) + b(4)))./(-a(1) + a(4)));expr130 = angle(1 + (i.*(b(1) + b(4)))./(-a(1) + a(4)));expr131 = angle(1 + (i.*(-b(2) + b(4)))./(-a(2) + a(4)));expr132 = angle(1 + (i.*(b(2) + b(4)))./(-a(2) + a(4)));expr133 = angle(1 + (i.*(-b(3) + b(4)))./(-a(3) + a(4)));expr134 = angle(1 + (i.*(b(3) + b(4)))./(-a(3) + a(4)));expr135 = -angle(1 + (i.*(b(4) - d(1)))./(a(4) - c(1)));expr136 = -angle(1 + (i.*(b(4) + d(1)))./(a(4) - c(1)));expr137 = -angle(1 + (i.*(b(4) - d(2)))./(a(4) - c(2)));expr138 = -angle(1 + (i.*(b(4) + d(2)))./(a(4) - c(2)));expr139 = -angle(1 + (i.*(b(4) - d(3)))./(a(4) - c(3)));expr140 = -angle(1 + (i.*(b(4) + d(3)))./(a(4) - c(3)));expr141 = -angle(1 + (i.*(b(4) - d(4)))./(a(4) - c(4)));expr142 = -angle(1 + (i.*(b(4) + d(4)))./(a(4) - c(4)));expr143 = a(4);expr144 = exp(x.*a(4));expr145 = cos(x.*b(4));expr146 = sin(x.*b(4));expr147 = (expr114.*expr115.*expr116.*expr117.*expr118.*expr119).^expr11;expr148 = (expr120.*expr121.*expr122.*expr123.*expr124.*expr125.*expr126.*expr127).^expr24;expr149 = expr128.*cos(expr129 + expr130 + expr131 + expr132 + expr133 + expr134 + expr135 + expr136 + expr137 + expr138 + expr139 + expr140 + expr141 + expr142) + expr143.*sin(expr129 + expr130 + expr131 + expr132 + expr133 + expr134 + expr135 + expr136 + expr137 + expr138 + expr139 + expr140 + expr141 + expr142);expr150 = sin(expr129 + expr130 + expr131 + expr132 + expr133 + expr134 + expr135 + expr136 + expr137 + expr138 + expr139 + expr140 + expr141 + expr142);expr151 = (expr128.*cos(expr129 + expr130 + expr131 + expr132 + expr133 + expr134 + expr135 + expr136 + expr137 + expr138 + expr139 + expr140 + expr141 + expr142) + expr143.*sin(expr129 + expr130 + expr131 + expr132 + expr133 + expr134 + expr135 + expr136 + expr137 + expr138 + expr139 + expr140 + expr141 + expr142)).^(-1);expr152 = (expr16.*expr17.*expr18.*expr19.*expr20.*expr21.*expr22.*expr23).^expr24;expr153 = (expr10.*expr5.*expr6.*expr7.*expr8.*expr9).^expr11;expr154 = expr25.*cos(expr26 + expr27 + expr28 + expr29 + expr30 + expr31 + expr32 + expr33 + expr34 + expr35 + expr36 + expr37 + expr38 + expr39) + expr40.*sin(expr26 + expr27 + expr28 + expr29 + expr30 + expr31 + expr32 + expr33 + expr34 + expr35 + expr36 + expr37 + expr38 + expr39);expr155 = sin(expr26 + expr27 + expr28 + expr29 + expr30 + expr31 + expr32 + expr33 + expr34 + expr35 + expr36 + expr37 + expr38 + expr39);expr156 = (expr25.*cos(expr26 + expr27 + expr28 + expr29 + expr30 + expr31 + expr32 + expr33 + expr34 + expr35 + expr36 + expr37 + expr38 + expr39) + expr40.*sin(expr26 + expr27 + expr28 + expr29 + expr30 + expr31 + expr32 + expr33 + expr34 + expr35 + expr36 + expr37 + expr38 + expr39)).^(-1);expr157 = (expr46.*expr47.*expr48.*expr49.*expr50.*expr51).^expr11;expr158 = (expr52.*expr53.*expr54.*expr55.*expr56.*expr57.*expr58.*expr59).^expr24;expr159 = expr60.*cos(expr61 + expr62 + expr63 + expr64 + expr65 + expr66 + expr67 + expr68 + expr69 + expr70 + expr71 + expr72 + expr73 + expr74) + expr75.*sin(expr61 + expr62 + expr63 + expr64 + expr65 + expr66 + expr67 + expr68 + expr69 + expr70 + expr71 + expr72 + expr73 + expr74);expr160 = sin(expr61 + expr62 + expr63 + expr64 + expr65 + expr66 + expr67 + expr68 + expr69 + expr70 + expr71 + expr72 + expr73 + expr74);expr161 = (expr60.*cos(expr61 + expr62 + expr63 + expr64 + expr65 + expr66 + expr67 + expr68 + expr69 + expr70 + expr71 + expr72 + expr73 + expr74) + expr75.*sin(expr61 + expr62 + expr63 + expr64 + expr65 + expr66 + expr67 + expr68 + expr69 + expr70 + expr71 + expr72 + expr73 + expr74)).^(-1);expr162 = (expr80.*expr81.*expr82.*expr83.*expr84.*expr85).^expr11;expr163 = (expr86.*expr87.*expr88.*expr89.*expr90.*expr91.*expr92.*expr93).^expr24;expr164 = expr94.*cos(expr100 + expr101 + expr102 + expr103 + expr104 + expr105 + expr106 + expr107 + expr108 + expr95 + expr96 + expr97 + expr98 + expr99) + expr109.*sin(expr100 + expr101 + expr102 + expr103 + expr104 + expr105 + expr106 + expr107 + expr108 + expr95 + expr96 + expr97 + expr98 + expr99);expr165 = sin(expr100 + expr101 + expr102 + expr103 + expr104 + expr105 + expr106 + expr107 + expr108 + expr95 + expr96 + expr97 + expr98 + expr99);expr166 = (expr94.*cos(expr100 + expr101 + expr102 + expr103 + expr104 + expr105 + expr106 + expr107 + expr108 + expr95 + expr96 + expr97 + expr98 + expr99) + expr109.*sin(expr100 + expr101 + expr102 + expr103 + expr104 + expr105 + expr106 + expr107 + expr108 + expr95 + expr96 + expr97 + expr98 + expr99)).^(-1);expr167 = expr113.*expr12.*expr13.*expr14.*expr147.*expr148.*expr149.*expr15.*expr2.*expr3.*(1 - expr144.*(expr145 - expr113.*expr146.*(expr143 - expr150.*expr151.*expr4))).*expr43;expr168 = expr1.*expr12.*expr13.*expr14.*expr15.*expr152.*expr153.*expr154.*expr2.*expr3.*expr4.*(1 - expr41.*(expr42 - expr1.*(expr40 - expr155.*expr156.*expr43).*expr44));expr169 = expr12.*expr13.*expr14.*expr15.*expr157.*expr158.*expr159.*expr3.*expr4.*expr43.*expr45.*(1 - expr76.*(expr77 - expr45.*(-(expr160.*expr161.*expr2) + expr75).*expr78));expr170 = expr12.*expr13.*expr14.*expr15.*expr162.*expr163.*expr164.*expr2.*expr4.*expr43.*expr79.*(1 - expr110.*(expr111 - expr112.*(expr109 - expr165.*expr166.*expr3).*expr79));expr171 = expr167 + expr168 + expr169 + expr170;f = M1.*expr171;end