clear all
close all
k=input('取樣點數=2^k,k=[14]');
if isempty(k),
  k=14;
end
N=2^k;  %采樣點數
L=64;  %每碼元采樣點數
M=N/L;  %碼元數
Rb=2;  %碼元速率
Ts=1/Rb;  %碼元間隔
dt=Ts/L;  %時域取樣間隔
df=1/(N*dt);  %頻域取樣間隔
T=N*dt;  %時域截短時間
Bs=N*df/2;  %頻域系統帶寬
t=[-T/2+dt/2:dt:T/2];  %時域橫坐標
f=[-Bs+df/2:df:Bs];    %頻域橫坐標 
%定義全零距陣，將來作為求功率的累加器
Ep=zeros(size(f));    
for ii=1:20
  n=1;  %賦初值
  while mod(n,2)~=0,  %保證1的個數為偶數個
    a=round(rand(1,M));
    n=length(find(a==1));
  end
  bb=find(a==1);  %找出a中1的位置
  %找出第奇數個1的位置
  cc=bb(1:2:length(bb));
  %將第奇數個位置的1變為-1  
  a(cc)=-ones(size(cc));  
  ami=zeros(L,M);
for loop=1:L
  ami(loop,:)=a;  %ami第loop行元素
end
ami=reshape(ami,1,N);
AMI=t2f(ami,dt);
P=AMI.*conj(AMI)/T;
Ep=(Ep*(ii-1)+P)/ii;
end
figure(1)
set(1,'position',[10,50,300,200])
aa=30+10*log10(Ep+eps);
plot(f,aa,'r')
grid on
axis([-5,5,-50,50])
xlabel('f(MHz)')
ylabel('Ps(f)(dBm/MHz)')
figure(2)
set(2,'position',[350,50,300,200])
plot(t,ami,'b')
grid on
axis([-5,5,-1.5,1.5])
xlabel('t(us)')
ylabel('s(t)(v)');