c10_mcqpskrun.m

QPSK调制下的蒙特-卡罗仿真主函数源程序

%file:c10_MCQPSKrun.m
function BER_MC=MCQPSKrun(N,Eb,N0,ChanAtt,TimingBias,Timingjitter,PhaseJitter)
fs=1e+6;                                                %sampling rate
SymRate=1e+5;                                           %symbol rate
Ts=1/fs;                                                %sampling period
TSym=1/SymRate;                                         %symbol period
SymToSend=N;                                            %symbol to be transmitted
ChanBW=4.99e+5;                                         %bandwidth of channel
MeanCarrierPhaseError=PhaseBias;                        %mean of carrier phase
StdCarrierPhaseError=PhaseJitter;                       %stdev of phase error
MeanSymbolSyncError=TimingBias;                         %mean of symbol sync error
StdSymbolSyncError=TimingJitter;                        %stdev of symbol sync error
ChanGain=10^(-ChanAtt/20);                              %channel gain(linear units)
TxBitClock=Ts/2;                                        %transmitter bit clock
RxBitClock=Ts/2;                                        %receive bit clock
%
%Standard deviation of noise and signal amplitude at receiver input.
%
RxNoiseStd=sqrt((10^((N0-30)/10)*(fs/2));               %stdeve of noise
TxSigAmp=sqrt(10^((Eb-30)/10)*SymRate);                 %signal amplitude
%
%Allocate some memory for probes.
%
SampPerSym=fs/SymRate;
probe1=zero((SymToSend+1)*SampPerSym,1);
probe1counter=1;
probe2=zero((SymToSend+1)*SampPerSym,1);
probe2counter=1;
%
%Counters to keep track of how many symbols have have been sent.
%
TxSymSent=1;
RxSymDemod=0;
%
%Buffers that contain the transmitted and received data.
%
[unused,SourceBitsI]=random_binary(SymToSend,1);
[unused,SourceBitsQ]=random_binary(SymToSend,1);
%
%Differentially encode the transmitted data.
%
TxBitsI=SourceBitsI*0;                                  %set first I bit
TxBitsQ=SourceBitsQ*0;                                  %set first Q bit
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
%
%Make a complex data stream of the I and Q bits.
%
TxBits=((TxBitsI*2)-1)+(sqrt(-1)*((TxBitsQ*2)-1));
%
RxIntegrator=0;
TxBitClock=2*TSym;
%
%Design the channel filter,and create the filter state array.
%
[b,a]=butter(2,ChanBW/(fs/2));
b=[1];a=[1];                                            %filter bypassed
[junk,FilterState]=filter(b,a,0);
%
%begin simulation loop.
%
while TxSymSent<SymToSend
    %
    %updata the transmitter's clock,and see if it is time to get new bits
    %
    TxBitClock=TxBitClock+Ts;
    if TxBitClock>TSym
        %
        %time to get new bits
        %
        TxSymSent=TxSymSent+1;
        %
        %we don't want the clock to increase to infinity,so subtract off
        %an integer number of Tb seconds.
        %
        TxBitClock=mod(TxBitClock,TSym);
        %
        %get the new bit,and scale it up appropriately.
        %
        TxOutput=TxBits(TxSymSent)*TxSigAmp;
    end
    %
    %pass the transmitted signal through the channel filter.
    %
    [Rx,FilterState]=filter(b,a,TxOutput,FilterState);
    %
    %add white gaussian noise to the signal.
    %
    Rx=(ChanGain*Rx)+(RxNoiseStd*(randn(1,1)+sqrt(-1)*randn(1,1)));
    %
    %phase rotation due to receiver carrier synchronization error.
    %
    PhaseRotation=exp(sqrt(-1)*2*pi*...
        (MeanCarrierPhaseError+(randn(1,1)*StdCarrierPhaseError))/360);
    Rx=Rx*PhaseRotation;
    probe1(probe1counter)=Rx; probe1counter=probe1counter+1;
    %
    %update the integrate and dump filter at the receiver.
    %
    RxIntegrator=RxIntegrator+Rx;
    probe2(probe2counter)=RxIntegrator;
    probe2counter=probe2counter+1;
    %
    %update the receiver clock,to see if it is time to sample and dump the
    %integrator.
    %
    RxBitClock=RxBitClock+Ts;
    xTSym=TSym*(1+MeanSymbolSyncError+(StdSymbolSyncError*randn(1,1)));
    if RxBitClock>RxTSym                                %time to demodulate symbol   
        RxSymDemod=RxSymDemod+1;
        RxBitsI(RxSymDemod)=round(sign(real(RxIntegrator))+1)/2;
        RxBitsQ(RxSymDemod)=round(sign(imag(RxIntegrator))+1)/2;
        RxBitClock=RxBitClock-TSym;                     %reset receive clock
        RxIntegrator=0;                                 %reset integrator
    end
end
%
%differential decoder.
%
SinkBitsI=SourceBitsI*0;                                %set first I sink bit
SinkBItsQ=SourceBitsQ*0;                                %set first Q sink bit
%
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;
%
%look for best time delay between input and output for 100 bits.
%
[C,Lags]=vxcorr(SourceBitsI(10,110),SinkBitsI(10,110));
[MaxC,LocMaxC]=max(C);
BestLag=Lags(LocMaxC);
%
%Adjust time delay to match best lag
%
if BestLag>0
    SourceBitsI=SourceBitsI(BestLag+1:length(SourceBitsI));
    SourceBitsQ=SourceBitsQ(BestLag+1:length(SourceBitsQ));
else if BestLag<0
    SinkBitsI=SinkBitsI(-BestLag+1:length(SinkBitsI));
    SinkBitsQ=SinkBitsQ(-BestLag+1:length(SinkBitsQ));
end
%
%make all arrays the same length.
%
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);
%
%find the number of errors and the BER.
%
Errors=sum(SourceBitsI~=SinkBitsI)+sum(SourceBitsQ~=SinkBitsQ);
BER_MC=Errors/2*length(SourceBItsI));
%end of function file