TxBitClock=TxBitClock+Ts; if TxBitClock>TSym
TxSymSent=TxSymSent+1;
TxBitClock=mod(TxBitClock,TSym);
TxOutput=TxBits(TxSymSent)*TxSigAmp; end
%信号经过信道滤波器
[Rx,FilterState]=filter(b,a,TxOutput,FilterState); %加高斯白噪声
Rx=(ChanGain*Rx)+(RxNoiseStd*(randn(1,1)+sqrt(-1)*randn(1,1))); %基于接收机载波同步误差的相位旋转
PhaseRotation=exp(sqrt(-1)*2*pi*(MeanCarrierPhaseError+(randn(1,1)*StdCarrierPhaseError))/360); Rx=Rx*PhaseRotation; probe1(probe1counter)=Rx;
probe1counter= probe1counter+1; %更新接收机积分清除器 RxIntegrator=RxIntegrator+Rx;
probe2(probe2counter)=RxIntegrator; probe2counter= probe2counter+1;
%更新接收机时钟,判断是不是适合采样 RxBitClock=RxBitClock+Ts;
RxTSym=TSym*(1+MeanSymbolSyncError+(StdSymbolSyncError*randn(1,1))); if RxBitClock>RxTSym %解调信号 RxSymDemod=RxSymDemod+1;
RxBitsI(RxSymDemod)=round(sign(real(RxIntegrator))+1)/2; RxBitsQ(RxSymDemod)=round(sign(imag(RxIntegrator))+1)/2; RxBitClock=RxBitClock-TSym; RxIntegrator=0; end end
% 差分解码
SinkBitsI=SourceBitsI*0; SinkBitsQ=SourceBitsQ*0;
for k=2:RxSymDemod
SinkBitsI(k)=or(and(not(xor(RxBitsI(k),RxBitsQ(k))),xor(RxBitsI(k),RxBitsI(k-1))),and(xor(RxBitsI(k),RxBitsQ(k)),xor(RxBitsQ(k),RxBitsQ(k-1))));
SinkBitsQ(k)=or(and(not(xor(RxBitsI(k),RxBitsQ(k))),xor(RxBitsQ(k),RxBitsQ(k-1))),and(xor(RxBitsI(k),RxBitsQ(k)),xor(RxBitsI(k),RxBitsI(k-1))));
end
%在输入和输出100字节中寻找最佳时延
[C,Lags]=vxcorr(SourceBitsI(10:110),SinkBitsI(10:110)); [MaxC,LocMaxC]=max(C); BestLag=Lags(LocMaxC); % 调整时延 if BestLag>0
SourceBitsI=SourceBitsI(BestLag+1:length(SourceBitsI)); SourceBitsQ=SourceBitsQ(BestLag+1:length(SourceBitsQ)); elseif BestLag<0
SinkBitsI=SinkBitsI(-BestLag+1:length(SinkBitsI)); SinkBitsQ=SinkBitsQ(-BestLag+1:length(SinkBitsQ)); end
% 将序列调整成相同长度
TotalBits=min(length(SourceBitsI),length(SinkBitsI)); TotalBits=TotalBits-20;
SourceBitsI=SourceBitsI(10:TotalBits); SourceBitsQ=SourceBitsQ(10:TotalBits); SinkBitsI=SinkBitsI(10:TotalBits); SinkBitsQ=SinkBitsQ(10:TotalBits);
Errors=sum(SourceBitsI ~= SinkBitsI)+sum(SourceBitsQ ~= SinkBitsQ); BER_MC=Errors/(2*length(SourceBitsI));
┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ 装 ┊ ┊ ┊ ┊ ┊ 订 ┊ ┊ ┊ ┊ ┊ 线 ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊
return
function MCQPSKBER
Eb=22:0.5:26; N0=-50;
ChannelAttenuation=70;
EbN0dB=(Eb-ChannelAttenuation)-N0; EbN0=10.^(EbN0dB./10); BER_T=0.5*erfc(sqrt(EbN0)); N=round(100./BER_T); BER_MC=zeros(size(Eb)); for k=1:length(Eb)
BER_MC(k)=MCQPSKrun(N(k),Eb(k),N0,ChannelAttenuation,0,0,0,0); disp(['simulation',num2str(k*100/length(Eb)),'%complate']); end
figure(1);
semilogy(EbN0dB,BER_MC,'o',EbN0dB,2*BER_T,'-'); xlabel('Eb/N0(dB)'); ylabel('bit error rate');
legend('MC BER Estimate','Theoretical BER');grid; return
代码3
function [peideal,pesystem]=qpsk_berest(xx,yy,ebn0db,eb,tb,nbw)
%xx为输入序列,yy为输出序列,ebn0db为信噪比,eb为单个信号能力,tb信号周期,nbw噪声带宽 [n1 n2]=size(xx); nx=n1*n2;
[n3 n4]=size(yy); ny=n3*n4;
[n5 n6]=size(ebn0db); neb=n5*n6;
%接收机带宽设定为rs/2 nbwideal=1/(2*tb*2); for m=1:neb
peideal(m)=0.0; % 初始化 pesystem(m)=0.0; % 初始化
%计算n0和噪声方差
string1=['Eb/N0=',num2str(ebn0db(m))]; disp(string1)
ebn0(m)=10^(ebn0db(m)/10);
n0=eb/ebn0(m); % 噪声功率 sigma=sqrt(n0*nbw*2); % 方差 sigma1=sqrt(n0*nbwideal*2); % 理想方差
b=sqrt(2*eb/tb)/sqrt(sum(abs(xx).^2)/nx); for n=1:nx
theta=angle(xx(n)); if(theta<0)
theta=theta+2*pi; end
%接收信号的相位旋转
xxx(n)=b*xx(n)*exp(-i*(theta-(pi/4))); yyy(n)=yy(n)*exp(-i*(theta-(pi/4))); d1=real(xxx(n)); d2=imag(xxx(n)); d3=real(yyy(n)); d4=imag(yyy(n));
┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ 装 ┊ ┊ ┊ ┊ ┊ 订 ┊ ┊ ┊ ┊ ┊ 线 ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊
pe1=q(d1/sigma1)+q(d2/sigma1); pe2=q(d3/sigma)+q(d4/sigma); peideal(m)=peideal(m)+pe1; pesystem(m)=pesystem(m)+pe2; end end
peideal=(1/2)*peideal./nx; pesystem=(1/2)*pesystem./nx; return
function QPSKThreeray(p0,p1,p2,delay)
NN=256; %发送的信号数 tb=0.5; %信号周期
fs=16; %每信号的采样数 ebn0db=[1:2:14]; %信噪比矢量 %产生QPSK信号
x=random_binary(NN,fs)+i*random_binary(NN,fs);
delay0=0; delay1=0; delay2=delay;
%设置瑞利信道
gain1=sqrt(p1)*abs(randn(1,NN)+i*randn(1,NN)); gain2=sqrt(p2)*abs(randn(1,NN)+i*randn(1,NN));
for k=1:NN for kk=1:fs
index=(k-1)*fs+kk;
ggain1(1,index)=gain1(1,k); ggain2(1,index)=gain2(1,k); end end y1=x;
for k=1:delay2
y2(1,k)=y1(1,k)*sqrt(p0); end
for k=(delay2+1):(NN*fs)
y2(1,k)=y1(1,k)*sqrt(p0)+y1(1,k-delay1)*ggain1(1,k)+y1(1,k-delay2)*ggain2(1,k); end
%匹配滤波器 b=-ones(1,fs); b=b/fs; a=1;
y=filter(b,a,y2);
%应用半解析方法,进行误码率估计。 [cor lags]=vxcorr(x,y); cmax=max(max(abs(cor))); nmax=find(abs(cor)==cmax); timelag=lags(nmax); corrmag=cmax;
theta=angle(cor(nmax)); y=y*exp(-i*theta);
hh=impz(b,a); ts=1/fs;
nbw=(fs/2)*sum(hh.^2);
index=(10*fs+8:fs:(NN-10)*fs+8); xx=x(index);
yy=y(index-timelag+1); [n1 n2]=size(y2);