%例二 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,
   a=round(rand(1,M));
   n=length(find(a==1));
 end
 bb=find(a==1);    %bb是a中1的位置
 cc=bb(1:2:length(bb)); 
    %cc是a中奇序數1的位置
 b=a;              
 b(cc)=-ones(size(cc)); %所有位于cc的1變成-1

 ami=zeros(L,M);
 for loop=1:L,ami(loop,:)=b; end
 ami=reshape(ami,1,N);
 
 AMI=t2f(ami);       %AMI信號的傅氏變換
 P=AMI.*conj(AMI)/T;        %AMI信號的功率譜
 EP=(EP*(ii-1)+P)/ii;       %累計平均

 figure(1)
 aa=30+10*log10(EP+eps);%加eps以避免除以零
 plot(f,aa,'r');
 grid
 axis([-5,+5,-50,50])
 xlabel('f  (MHz)')
 ylabel('Ps(f)   (dBm/MHz)')

 figure(2)
 plot(t,ami,'r')
 grid
 axis([-2,+2,-1.5,1.5])
 xlabel('t (us)')
 ylabel('s(t) (V)')

end