c22.m

来自「数值模拟实验」· M 代码 · 共 65 行

function D=c22(k)
W=zeros(1,2^k);
C=zeros(1,2^k/2);
S=zeros(1,2^k/2);
for l=1:2^k/2
    if l==1,
        C(l)=1;
        S(l)=0;
    else if l==2,
            C(l)=cos(2*pi/(2^k));
            S(l)=sin(2*pi/(2^k));
        else
            C(l)=C(2)*C(l-1)-S(2)*S(l-1);
            S(l)=S(2)*C(l-1)+C(2)*S(l-1);
        end
    end
    W(l)=C(l)+i*S(l);
    W(l+2^(k-1))=-W(l);
end
B_0=c21(k);  
n1=zeros(1,2^k);  
Order=zeros(1,k);
for n2=1:2^k
    n3=n2-1;
    for h=1:k   %change the order
        if n3>=0
           Order(h)=rem(n3,2);
           n3=(n3-Order(h))/2;
           n1(n2)=n1(n2)+2^(k-h)*Order(h);   %changed  
        end
    end
end
%%%%%
B=zeros(2^k,k+1);
for m=1:2^k
    B(m,1)=B_0(n1(m)+1);
end
for row=2:k+1
    rr=row-1;                  %
    for list=1:2^k 
        order1=n1(list);  %
        for tt=1:(k-rr+1)
            rea=rem(order1,2);
            order1=(order1-rea)/2;
        end
        order2=2^(k-rr)*rem((list-1),2^rr)+1;      %
        if  rea==1
            value=n1(list)-2^(k-rr);
            for n4=1:2^k
                if n1(n4)==value
                    n5=n4;
                end
            end
            B(rr*2^k+list)=B((rr-1)*2^k+list)*W(order2)+B((rr-1)*2^k+n5);
        else  value=n1(list)+2^(k-rr);
            for n4=1:2^k
                if n1(n4)==value
                    n5=n4;
                end
            end
            B(rr*2^k+list)=B((rr-1)*2^k+n5)*W(order2)+B((rr-1)*2^k+list);
        end
    end
end
D=B(:,k+1);