lxk5.m

本程序仿真一个从卫星到区域数据中心的下行链路（由一组6个载波构成）

a1=1.0;a2=1.0;a3=1.0;
a4=1.0;a5=0.0;a6=1.0;
ibo=-5;
R=1;tsymbol=1/R;
df=1.244;
beta=0.2;
nsamples=64;
nsymbols=256;
ncorr=64;
nstart=64;
ebn0db=0:1:8;
tb=tsymbol/2;fs=nsamples*R;ts=1/fs;
f1=(df/2);f2=f1+df;f3=f2+df;f4=-f1;f5=-f2;f6=-f3;
delay1=0;delay2=8;delay3=16;
delay4=32;delay5=40;delay6=50;
phase1=0.0;phase2=pi/64;phase3=pi/32;
phase4=pi/48;phase5=pi/16;phase6=pi/40;
omega1=2*pi*f1;omega2=2*pi*f2;omega3=2*pi*f3;
omega4=-omega1;omega5=-omega2;omega6=-omega3;
M=4;
[xin1]=mpsk_impulses(M,nsymbols,nsamples);
b=ones(1,nsamples);
x=filter(b,1,xin1);
nduration=8;
[imp1,impw]=sqrc_time(R,beta,nsamples,nduration);
txout1=filter(imp1,1,xin1);
[xin2]=mpsk_impulses(M,nsymbols,nsamples);
[xin3]=mpsk_impulses(M,nsymbols,nsamples);
[xin4]=mpsk_impulses(M,nsymbols,nsamples);
[xin5]=mpsk_impulses(M,nsymbols,nsamples);
[xin6]=mpsk_impulses(M,nsymbols,nsamples);
[xin2]=delayr1(xin2,delay2);
[xin3]=delayr1(xin3,delay3);
[xin4]=delayr1(xin4,delay4);
[xin5]=delayr1(xin5,delay5);
[xin6]=delayr1(xin6,delay6);
txout2=filter(imp1,1,xin2);
txout3=filter(imp1,1,xin3);
txout4=filter(imp1,1,xin4);
txout5=filter(imp1,1,xin5);
txout6=filter(imp1,1,xin6);
KNN=1:nsymbols*nsamples;
txout1=a1*txout1.*exp(i*omega1*KNN*ts)*exp(i*phase1);
txout2=a2*txout2.*exp(i*omega2*KNN*ts)*exp(i*phase2);
txout3=a3*txout3.*exp(i*omega3*KNN*ts)*exp(i*phase3);
txout4=a4*txout4.*exp(i*omega4*KNN*ts)*exp(i*phase4);
txout5=a5*txout5.*exp(i*omega5*KNN*ts)*exp(i*phase5);
txout6=a6*txout6.*exp(i*omega6*KNN*ts)*exp(i*phase6);
txout=txout1+txout2+txout3+txout4+txout5+txout6;
[bmi ami]=butter(4,4.5/(fs/2));
imuxout=filter(bmi,ami,txout);
run('twt_data');twtdata=data;
ampout=twt_model(imuxout,twtdata,ibo);
[bmo amo]=butter(4,4.5/(fs/2));
omuxout=filter(bmo,amo,ampout);
psdflag=1;
if psdflag==1
    nfft=(nsymbols/2)*nsamples;nmax=nsymbols*nsamples;
    temp(1:nfft)=imuxout(nfft+1:nmax);
    [logpsd,freq,ptotal,pmax]=log_psd(temp,4096,ts);
    figure;subplot(2,1,1);
    plot(freq,logpsd);axis([-10 10 -40 0]);
    title('PSD of the TWTA input');ylabel('PSD');grid;
    temp(1:nfft)=ampout(nfft+1:nmax);
    [logpsd,freq,ptotal,pmax]=log_psd(temp,4096,ts);
    subplot(2,1,2)
    plot(freq,logpsd);axis([-10 10 -40 0]);
    title('PSD of the TWTA output');ylabel('PSD');grid;
end
[n1 n2]=size(omuxout);nxx=n1*n2;
eb=0.5*tb*sum(sum(abs(omuxout).^2))/nxx;
eb=eb*(a1*a1)/(a1*a1+a2*a2+a3*a3+a4*a4+a5*a5+a6*a6)
ydemod=omuxout.*exp(-i*omega1*KNN*ts);
nduration=16;
[imp2,impw]=sqrc_freq_nosinc(R,beta,nsamples,nduration);
y=filter(imp2,1,ydemod);
hh=impz(imp2,1);ts=1/fs;
nbw=(fs/2)*sum(abs(hh).^2)
corlength=ncorr*nsamples;
[cor lags]=vxcorr(x(1:corlength),y(1:corlength));
cmax=max(max(abs(cor)));nmax=find(abs(cor)==cmax);
timelag=-lags(nmax);corrmag=cmax;theta=-angle(cor(nmax));
y=y*exp(-i*theta);
maxindex=nsymbols*nsamples-(3*nsamples/2)-timelag;
index=((nstart)*nsamples+(nsamples/2):nsamples:maxindex);
timelag1=timelag+(nsamples/2);
xx=x(index);
yy=y(index+timelag1+1);
[peideal,pesystem]=qpsk_berest(xx,yy,ebn0db,eb,tb,nbw);
figure;subplot(1,2,1)
yscale=1.5*max(real(yy));
plot(yy,'+')
xlabel('Direct Sample');ylabel('Quadrature Sample');grid;
axis([-yscale yscale -yscale yscale])
subplot(1,2,2)
semilogy(ebn0db,peideal,'-.',ebn0db,pesystem);grid;
xlabel('E_b/N_0');ylabel('Bit Error Rate')
legend('Ideal','System')