part.m

生成DCT矩阵,只要简单输入矩阵的阶数,就可以返回一个用语DCT变换的矩阵

g=ones(40,1);
ff=optimset;
ff.Display='iter';
ff.MaxFunEvals=1e10;
ff.MaxIter=1e15;
ff.TolX=1e-20;
syms w z;
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)');
x41=sym('x(41)');
x42=sym('x(42)');
x43=sym('x(43)');
x44=sym('x(44)');
x45=sym('x(45)');
x46=sym('x(46)');%程序前期设置



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];
b=inv([1 0 0;x4 1 0;x5 x6 1]*diag([x1 x2 x3])*[1 x7 x8;0 1 x9;0 0 1]);
k=inv([1 0 0;x16 1 0;x17 x18 1]*diag([x13 x14 x15])*[1 x19 x20;0 1 x21;0 0 1]);

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;

invu0=([eye(3) 0]+[x41;x42;x43]*[[-x10 -x11 -x12] eye(3)])*diag(b,1);
invv0=([eye(3) 0]+[x44;x45;x46]*[[-x22 -x23 -x24] eye(3)])*diag(k,1);

invE0=1/2^(1/2)*[eye(3),zeros(3);zeros(3),[0 0 1;0 1 0;1 0 0]]*[eye(3),eye(3);eye(3),-eye(3)]*diag([invu0 invv0]);
invG1=1/2*[eye(4),eye(4);eye(4),-eye(4)]*[eye(4),zeros(4);zeros(4),z*eye(4)]*[eye(4),eye(4);eye(4),-eye(4)]*[inv(u1),zeros(4);zeros(4),eye(4)];
R=z^(-1)*invE0*invG1;

h=subs(E,z,z^6)*[1;z^(-1);z^(-2);z^(-3);z^(-4);z^(-5)];
H1=h(1)
H2=h(2)
H3=h(3)
H4=h(4)
H5=h(5)
H6=h(6)
H7=h(7)
H8=h(8)