%例二 AMI碼的功率譜

%本程序中時間單位是微秒
%頻率單位為MHz
%碼元速率單位是Mb/s

global dt df t f N
close all

N=2^13;   %采樣點數
L=64;    %每碼元的采樣點數
M=N/L     %碼元數
Rb=2;
Ts=0.5		%碼元寬度是0.5us
dt=Ts/L;        
df=1/(N*dt)     %MHz
T=N*dt          %截短時間
Bs=N*df/2       %系統帶寬

t=[-T/2+dt/2:dt:T/2]; %時域橫坐標
f=[-Bs+df/2:df:Bs];   %頻域橫坐標

%figure(1)
 %set(1,'Position',[10,50,300,200])
                %設定窗口位置及大小 
%figure(2)
 %set(2,'Position',[400,50,300,200])
                %設定窗口位置及大小 

EP=zeros(size(f));
for ii=1:50
 n=1;  %每次有偶數個1
 while rem(n,2)~=0,     %求余數，n為奇數
   a=round(rand(1,M));  %M個隨機的01碼
   n=length(find(a==1));    %1的個數
 end
 bb=find(a==1);    %bb是a中1的位置，例如bb=1 2 3 4 9 10
 %cc=bb(1:2:length(bb)); 
    %cc是a中奇序數1的位置，若bb如上，則cc=1 3 9
 b=a;     %b，a相當于一個數組
 %b(cc)=-ones(size(cc)); %所有位于cc的1變成-1

 ami=zeros(L,M);    %L行M列的零矩陣
 for loop=1:L,ami(loop,:)=b; end  %每一行都賦為a的值，M列即M個碼元，L行即每碼元采樣L點
 for loop=L/2+1:L,ami(loop,:)=zeros(1,M); end  %所有碼元后半周期歸零
 ami=reshape(ami,1,N);  %將ami轉成1行N列，以便作為橫坐標
 
 AMI=t2f(ami);       %AMI信號的傅氏變換
 P=AMI.*conj(AMI)/T;        %AMI信號的功率譜；conj是取復共軛
 EP=(EP*(ii-1)+P)/ii;       %累計平均
                            %length（EP）/M=？
 
 %figure(1)
 

 %figure(2)
 subplot(2,1,1)
 plot(t,ami,'r')
 grid
 axis([-3,+3,-1.5,1.5])
 xlabel('t (us)')
 ylabel('s(t) (V)')
 title('單極性歸零碼碼型（Ts=0.5us）')
 
 subplot(2,1,2)
 aa=30+10*log10(EP+eps);%加eps以避免除以零,如果把EP換成P，則圖形有毛刺；隨機碼越長，圖形越平滑
 plot(f,aa,'r');
 grid
 axis([-15,+15,-50,50])
 xlabel('f  (MHz)')
 ylabel('Ps(f)   (dBm/MHz)')
 title('單極性歸零碼功率譜')

end