qpsk_awgn_ber.m.txt

AWGN多径信道中的QPSK系统进行BER性能仿真

%File: c14_threeray.m
%%%%%%%%%%%%%%%%
%Default parameters

NN = 256;  %number of symbols
tb = 0.5;  %bit time
fs = 16;   %samples/symbol
ebn0db = [1:2:14];  %Eb/N0 vector

%%%%%%%%%%%%%%%%%%%%%%%%
%Establish QPSK signals
%%%%%%%%%%%%%%%%%%%%%%%%%%

x=random_binary(NN,fs)+i*random_binary(NN,fs);   %QPSK signal

%%%%%%%%%%%%%%%%%%%%%%%
%Input power and delays
%%%%%%%%%%%%%%%%%%%%%%%%%
p0 = input('Enter p0>');
p1 = input('Enter p1>');
p2 = input('Enter p2>');
delay = input('Enter tau>');
delay0 = 0;
delay1 = 0;
delay2 = delay;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set up the Complex Gaussian (Rayleigh) gains'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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

%%%%%%%%%%%%%%%%%%%%%%%
% Matched filter
%%%%%%%%%%%%%%%%%%%%%%%
b = -ones(1,fs);
b = b/fs;
a = 1;
y = filter(b,a,y2);

%%%%%%%%%%%%%%%%%%%%%%%%%%
%End of simulation
%%%%%%%%%%%%%%%%%%%%%%%%
%Use the semi analytic BER estimator.The following sets
%up the semi analytic estimator. Find the maximum magnitude 
%of the cross correlation and the corresponding lag.
%%%%%%%%%%%%%%%%%%%%%%%%
[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);  %derotate

%%%%%%%%%%%%
%Noise BW calibration
%%%%%%%%%%%%%
hh = impz(b,a);
ts = 1/16;
nbw = (fs/2)*sum(hh.^2);

%%%%%%%%%%%%%%
%Delay the input,and do BER estimation on the last 128bits.
%Use middle sample. Make sure the index does not exceed number
%of input points.Eb should be computed at the receiver input.
%%%%%%%%%%%%%
index = (10*fs+8:fs:(NN-10)*fs+8);
xx = x(index);
yy = y(index-timelag+1);
[n1 n2] = size(y2);
ny2=n1*n2;
eb = tb*sum(sum(abs(y2).^2))/ny2;
eb = eb/2;
[peideal,pesystem] = qpsk_berest(xx,yy,ebn0db,eb,tb,nbw);
figure;
semilogy(ebn0db,peideal,'b*-',ebn0db,pesystem,'r+-');
xlabel('E_b/N-0(dB)');
ylabel('Probability of Error');
grid;
axis([a 14 10^(-10) 1]);

%End of script file.


%file: random_binary.m
function [x,bits] = random_binary(nbits,nsamples)
%This function genrates a random binary waveform of length nbits
%sample at a rate of nsamples/bit.
x = zeros(1,nbits*nsamples);
bits = round(rand(1,nbits));
for m=1:nbits
    for n=1:nsamples
       index = (m-1)*nsamples + n;
       x(1,index) = (-1)^bits(m);
      end
   end
%End of function file.


%File: vxcorr.m
function [c,lags] = vxcorr(a,b)
%This function calculates the unscaled cross-correlation of 2
%vectors of the same length. The output length(c) is length(a)+length(b)-1.
%It is a simplified function of xcorr function in matlabR12 using the defination:
%c(m) = E[a(n+m)*conj(b(n))] = E[a(n)*conj(b(n-m))]
%
a=a(:);     %convert a to column vector
b=b(:);    %convert b to column vector
M=length(a);   %same as length(b)
maxlag = M-1;   %maximum value of lag
lags = [-maxlag:maxlag]';    %vector of lags
A=fft(a,2^nextpow2(2*M-1));   %fft od A
B=fft(b,2^nextpow2(2*M-1));   %fft of B
c=ifft(A.*conj(B));   %crosscorrelation
%
%Move negative lags before positive lags.
%
c = [c(end-maxlag+1:end,1);c(1:maxlag+1,1)];
%
%Return row vector if a,b are row vectors.

[nr nc]=size(a);
if(nr>nc)
   c=c.';
   lags=lags.';
end
%End of function file.

%File: psk_berest.m
function [peideal,pesystem]=psk_berest(xx,yy,enn0db,eb,tb,nbw)
%ebn0db is an array of Eb/N0 values in db(specified at the receiver input);
%tb is the bit duration and nbw is the noise BW
%xx is the reference (ideal)input ;
%
nx = length(xx);
%
%For comparision purposes, set the noise BW of the ideal
%receiver (integrate and dump) to be equal to rs/2.
%
nbwideal = 1/(2*tb);    %noise bandwidth
for m=1:length(ebn0db)
    peideal(m)=0.0;  pesystem(m)=0.0   %initialize
    %
    %Find n0 and the variance of the noise.
    %
    ebn0(m) = 10^(ebn0db(m)/10);     %dB to linear
    n0 = eb/ebn0(m);      %noise power
    sigma = sqrt(n0*nbw*2);   %variance
    sigma1 = sqrt(n0*nbwideal*2);   %variance of ideal
    % 
    %Multiply the input constellation/signal by a scale factor so
    %that input constellation and the constellations/signal at the
    %input to receive filter have the same ave power
    % a = sqrt(2*eb/(2*tb))
    %
    b = sqrt(2*eb/tb)/sqrt(sum(abs(xx).^2)/nx);
    d1 = b*abs(xx);
    d3 = abs(yy);
    peideal(m) = sum (q(d1/sigma1));
    pesystem(m) = sum (q(d3/sigma));
  end
peideal = peideal/nx;
pesystem = pesystem/nx;
%End of function file.

%File: q.m
function out=q(x)
out = 0.5*erfc(x/sqrt(2));
%End of function file.