N=127;   %周期為N=127
a=2;    % 序列的幅值
ndp=254*2; %輸出的點數，如果ndp=254說明兩個周期的m序列
for(i=1:7)
    SH(i)=2*a;
end
for(i=1:ndp)
    s=SH(7)+SH(4);
    if(s>3*a)
        s=0.0;
    end
    for(j=2:7)
        jm=8-j;
        SH(jm+1)=SH(jm);
    end
    SH(1)=s;
    U(i)=s-a;          %U為輸出的m序列
end


%下面是逆重復m序列  《動、靜態數學模型的實用建模方法》P363---P369
for (i=1:2:ndp)
     U2(i)=-(U(i));
 end
 for (i=2:2:ndp)
     U2(i)=U(i);
 end
 figure(1);
 subplot(2,1,1);
 plot(U);
 subplot(2,1,2);
 plot(U2);

%m序列的自相關函數
for(i=0:ndp)
    he=0.0;
    for(j=1:ndp-i)
        he=he+U(j)*U(j+i);
    end
    r(i+1)=he/(ndp-i);   %自相關函數的第二種定義為r(i+1)=he/ndp；
end
%%逆重復m序列的自相關函數
for(i=0:ndp)
    he=0.0;
    for(j=1:ndp-i)
        he=he+U2(j)*U2(j+i);
    end
    r2(i+1)=he/(ndp-i);%自相關函數的第二種定義為r2(i+1)=he/ndp；
end
figure(2);
subplot(2,1,1);
plot(-254:254,r);
subplot(2,1,2);
plot(-254:254,r2);