analysis1.m

用于多通道滤波器组实现最优化设计,主要应用于图象处理

g=ones(40,1);
ff=optimset;
ff.Display='iter';
ff.MaxFunEvals=1e10;
ff.MaxIter=1e15;
ff.TolX=1e-20;

x1=sym('x(1)');
x2=sym('x(2)');
x3=sym('x(3)');
x4=sym('x(4)');
x5=sym('x(5)');
x6=sym('x(6)');
x7=sym('x(7)');
x8=sym('x(8)');
x9=sym('x(9)');
x10=sym('x(10)');
x11=sym('x(11)');
x12=sym('x(12)');

x13=sym('x(13)');
x14=sym('x(14)');
x15=sym('x(15)');
x16=sym('x(16)');
x17=sym('x(17)');
x18=sym('x(18)');
x19=sym('x(19)');
x20=sym('x(20)');
x21=sym('x(21)');
x22=sym('x(22)');
x23=sym('x(23)');
x24=sym('x(24)');
x25=sym('x(25)');
x26=sym('x(26)');
x27=sym('x(27)');
x28=sym('x(28)');
x29=sym('x(29)');
x30=sym('x(30)');
x31=sym('x(11)');
x32=sym('x(32)');
x33=sym('x(33)');
x34=sym('x(34)');
x35=sym('x(35)');
x36=sym('x(36)');
x37=sym('x(37)');
x38=sym('x(38)');
x39=sym('x(39)');
x40=sym('x(40)');         %程序前期设置



u0=[[1 0 0;x4 1 0;x5 x6 1]*diag([x1 x2 x3])*[1 x7 x8;0 1 x9;0 0 1];[x10 x11 x12]];
v0=[[1 0 0;x16 1 0;x17 x18 1]*diag([x13 x14 x15])*[1 x19 x20;0 1 x21;0 0 1];[x22 x23 x24]];
u1=[[1 0 0 0;x29 1 0 0;x30 x31 1 0;x32 x33 x34 1]*diag([x25 x26 x27 x28])*[1 x35 x36 x37;0 1 x38 x39;0 0 1 x40;0 0 0 1]]
J3=[0 0 1;0 1 0;1 0 0];


E0=1/2^(1/2)*[u0,zeros(4,3);zeros(4,3),v0]*[eye(3),eye(3);eye(3),-eye(3)]*[eye(3),zeros(3);zeros(3),[0 0 1;0 1 0;1 0 0]]
G1=1/2*[u1,zeros(4);zeros(4),eye(4)]*[eye(4),eye(4);eye(4),-eye(4)]*[eye(4),zeros(4);zeros(4),z^(-1)*eye(4)]*[eye(4),eye(4);eye(4),-eye(4)];
E=G1*E0;

E00=1/sqrt(8)*[u1*u0+u1*v0 u1*u0*J3-u1*v0*J3;u0+v0 u0*J3-v0*J3];                %E00表示E(z)的常数项
E01=1/sqrt(8)*[u1*u0-u1*v0 u1*u0*J3+u1*v0*J3;-u0+v0 -u0*J3-v0*J3];              %E01表示E(z)的z^(-1)项

a1=[0 1/24 1/4 1];
b1=[1 1 0 0];
P1=firpm(11,a1,b1);

a2=[0 2/24 3/24 5/24 8/24 1];
b2=[0 0 0.6 1 0 0];
P2=1.2*firpm(11,a2,b2);

a3=[0 5/24 6/24 8/24 12/24 1];
b3=b2;
P3=firpm(11,a3,b3);

a4=[0 4/24 9/24 11/24 14/24 1];
b4=b3;
P4=1.1*firpm(11,a4,b4);

a5=[0 6/24 14/24 16/24 18/24 1];
b5=[0 0 1 0.6 0 0];
P5=0.9*firpm(11,a5,b5,'differentiator');



a6=[0 8/24 17/24 19/24 21/24 1];
b6=[0 0 1 0.6 0 0];
P6=0.9*firpm(11,a6,b6,'differentiator');



a7=[0 12/24 20/24 21/24 22/24 1];
b7=[0 0 0.9 0.5 0 0];
P7=firls(11,[0 .1 .15 .2 .3 .4 .5 .6 .7 .8 .9 1],[0 0 0 0 0 0 0 0 1 1 0 0],'differentiator');


P8=firpm(11,[0 .4 .9 1],[0 0 1 1],'h');                %理想的8个滤波器


Ed0=[P1(1) P1(2) P1(3) P1(4) P1(5) P1(6);P2(1) P2(2) P2(3) P2(4) P2(5) P2(6);P3(1) P3(2) P3(3) P3(4) P3(5) P3(6);P4(1) P4(2) P4(3) P4(4) P4(5) P4(6);P5(1) P5(2) P5(3) P5(4) P5(5) P5(6);P6(1) P6(2) P6(3) P6(4) P6(5) P6(6);P7(1) P7(2) P7(3) P7(4) P7(5) P7(6);P8(1) P8(2) P8(3) P8(4) P8(5) P8(6)];
Ed1=[P1(7) P1(8) P1(9) P1(10) P1(11) P1(12);P2(7) P2(8) P2(9) P2(10) P2(11) P2(12);P3(7) P3(8) P3(9) P3(10) P3(11) P3(12);P4(7) P4(8) P4(9) P4(10) P4(11) P4(12);P5(7) P5(8) P5(9) P5(10) P5(11) P5(12);P6(7) P6(8) P6(9) P6(10) P6(11) P6(12);P7(7) P7(8) P7(9) P7(10) P7(11) P7(12);P8(7) P8(8) P8(9) P8(10) P8(11) P8(12)];
%提取8个理想滤波器组的系数

k1=[E00(1,:) E01(1,:)]
k2=[E00(2,:) E01(2,:)]
k3=[E00(3,:) E01(3,:)]
k4=[E00(4,:) E01(4,:)]
k5=[E00(5,:) E01(5,:)]
k6=[E00(6,:) E01(6,:)]
k7=[E00(7,:) E01(7,:)]
k8=[E00(8,:) E01(8,:)]                                  %提取8个设计的分析滤波器组的系数；

s=trace((k1-P1)*(k1-P1)')+trace((k2-P2)*(k2-P2)')+trace((k3-P3)*(k3-P3)')+trace((k4-P4)*(k4-P4)')+trace((k5-P5)*(k5-P5)')+2.3*trace((k6-P6)*(k6-P6)')+1.5*trace((k7-P7)*(k7-P7)')+15*trace((k8-P8)*(k8-P8)')
%优化分析滤波器组的目标函数