clear all;
close all;
%雙相碼波形及功率譜密度

L=16;%每個碼元的采樣點數
TS=1;%碼元周期
x=round(rand(1,1000));%產生0，1隨機序列
t=0:1/L:length(x);%定義對應時間序列
for i=1:length(x)%進行編碼
    if(x(i)==1)%若輸入信息為1
        for j=1:L/2;
            y(L/2*(2*i-2)+j)=1;%定義前半時間為1
            y(L/2*(2*i-1)+j)=-1;%定義后半時間為-1
        end
    else
        for j=1:L/2 %反之，輸入信息為0
            y(L/2*(2*i-2)+j)=-1;%定義前半時間為－1
             y(L/2*(2*i-1)+j)=1;%定義后半時間為1
        end;end;end
y=[y,x(i)];
figure(1)
subplot(221)
plot(t,y); 
axis([0 20 -1.5 1.5]);xlabel('雙相碼波形');
[f st1f]=T2F(t,y);
subplot(222)
plot(f,10*log10(abs(st1f+eps).^2/TS));grid
axis([-5 5 -40 20]);xlabel('f');ylabel('雙相碼功率譜密度(dB/Hz)');


%密勒碼波形及功率譜密度
t=0:1/L:length(x);
i=1; %由于第一碼元的編碼不定，直接給出
if(x(i)==1) %輸入為1，則
    for j=1:L/2;
        y2(L/2*(2*i-2)+j)=1;%前半時間為1
        y2(L/2*(2*i-1)+j)=-1;%后半時間為-1
    end
else                       %反之，輸入為0
    for j=1:L;
        y2(L*(i-1)+j)=1;%所有時間為1
    end
end
for i=2:length(x)%開始進行密勒編碼
    if(x(i)==1)  %輸入信息為1
        for j=1:L/2
            y2(L/2*(2*i-2)+j)=y2(L/2*(2*i-3)+L/4);%前半時間與前一碼元后半時間取值相同
            y2(L/2*(2*i-1)+j)=-y2(L/2*(2*i-2)+j);%后半時間與前半時間相反
        end
    else
        if(x(i-1)==1)%反之，若前一信息為1，輸入信息0
            for j=1:L
                y2(L*(i-1)+j)=y2(L/2*(2*i-3)+L/4);%所有時間與前一碼元后半時間值相同
            end
        else    %否則，前一信息為0
            for j=1:L
                y2(L*(i-1)+j)=-y2(L/2*(2*i-3)+L/4);%所有時間與前一碼元后半時間值相反
            end
        end
    end
end
y2=[y2 x(i)];
subplot(223)
plot(t,y2);
axis([0 20 -1.5 1.5]);
xlabel('Miller碼波形');
[f st2f]=T2F(t,y2);
subplot(224)
plot(f,10*log10(abs(st2f+eps).^2/TS));grid
axis([-5 5 -40 20]);xlabel('f');ylabel('Miller功率譜密度(dB/Hz)');