基于matlab的QPSK与BPSK信号性能比较仿真

┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ 装 ┊ ┊ ┊ ┊ ┊ 订 ┊ ┊ ┊ ┊ ┊ 线 ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊

function BER_MC=MCQPSKrun(N,Eb,N0,ChanAtt,TimingBias,TimingJitter,PhaseBias,PhaseJitter)

fs=1e+6; % 采样速率 SymRate=1e+5; % 信号速率 Ts=1/fs; % 采样周期 TSym=1/SymRate; % 信号周期 SymToSend=N; % 发送信号数 ChanBW=4.99e+5; % 带宽

MeanCarrierPhaseError=PhaseBias; % 载波相位均值 StdCarrierPhaseError=PhaseJitter; % 相位误差标准差

MeanSymbolSyncError=TimingBias; % 符号同步误差均值 StdSymbolSyncError=TimingJitter; % 符号同步误差标准差 ChanGain=10^(-ChanAtt/20); % 信道增益 TxBitClock=Ts/2; % 发送机时钟 RxBitClock=Ts/2; % 接收机时钟

%接收机输入端噪声标准差和信号幅度

RxNoiseStd=sqrt((10^((N0-30)/10))*(fs/2)); TxSigAmp=sqrt(10^((Eb-30)/10)*SymRate);

%分配缓存

SampPerSym=fs/SymRate;

probe1=zeros((SymToSend+1)*SampPerSym,1); probe1counter=1;

probe2=zeros((SymToSend+1)*SampPerSym,1); probe2counter=1;

%计已传输信号的个数

TxSymSent=1; RxSymDemod=0;

%发送和接收数据缓冲区

[unused,SourceBitsI]=random_binary(SymToSend,1); [unused,SourceBitsQ]=random_binary(SymToSend,1);

%差分编码

TxBitsI=SourceBitsI*0; TxBitsQ=SourceBitsQ*0; for k=2:length(TxBitsI)

TxBitsI(k)=or(and(not(xor(SourceBitsI(k),SourceBitsQ(k))),xor(SourceBitsI(k),TxBitsI(k-1))),and(xor(SourceBitsI(k),SourceBitsQ(k)),xor(SourceBitsQ(k),TxBitsQ(k-1))));

TxBitsQ(k)=or(and(not(xor(SourceBitsI(k),SourceBitsQ(k))),xor(SourceBitsQ(k),TxBitsQ(k-1))),and(xor(SourceBitsI(k),SourceBitsQ(k)),xor(SourceBitsI(k),TxBitsI(k-1))));

end

% 产生复信号

TxBits=((TxBitsI*2)-1)+(sqrt(-1)*((TxBitsQ*2)-1));

RxIntegrator=0; %初始化接收机积分器 TxBitClock=2*TSym; %初始化发送机

%设计信道滤波器，产生滤波器参数序列 [b,a]=butter(2,ChanBW/(fs/2)); b=[1]; a=[1];

[junk,FilterState]=filter(b,a,0);

% 开始仿真循环

┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ 装 ┊ ┊ ┊ ┊ ┊ 订 ┊ ┊ ┊ ┊ ┊ 线 ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊ ┊

while TxSymSent

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);