% Default parameters %QOSK
NN = 256; % number of symbols tb = 0.5; % bit time p0 = 1; % power
fs = 16; % samples/symbol ebn0db = [0:1:15]; % Eb/N0 vector
[b,a] = butter(5,1/16); % transmitter filter parameters %
% Establish QPSK signals %
x = random_binary(NN,fs)+1i*random_binary(NN,fs); % QPSK signal y1 = x; % save signal y2a = y1*sqrt(p0); % scale amplitude %
% Transmitter filter %
y2 = filter(b,a,y2a); % filtered signal %
% Matched filter %
b = ones(1,fs); b = b/fs; a = 1; % matched filter parameters y = filter(b,a,y2); % matched filter output %
% End of simulation %
% Use the semianalytic 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(abs(cor));
nmax = find(abs(cor)==cmax); timelag = lags(nmax); theta = angle(cor(nmax));
y = y*exp(-1i*theta); % derotate %
% Noise BW calibration %
hh = impz(b,a); % receiver impulse response
nbw = (fs/2)*sum(hh.^2); % noise bandwidth %
% Delay the input, and do BER estimation on the last 128 bits.
% 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); 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,'ro-'); grid; xlabel('E_b/N_0 (dB)'); ylabel('Bit Error Rate') legend('AWGN Reference','System Under Study') % End of script file.
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 % definition: 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 of 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.
function [peideal,pesystem] = qpsk_berest(xx,yy,ebn0db,eb,tb,nbw) % ebn0db is an array of Eb/No 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; yy is the distorted output; %
[n1 n2] = size(xx); nx = n1*n2; [n3 n4] = size(yy); ny = n3*n4;
[n5 n6] = size(ebn0db); neb = n5*n6; %
% For comparision purposes, set the noise BW of the ideal % receiver (integrate and dump) to be equal to rs/2. %
nbwideal = 1/(2*tb*2); for m=1:neb
peideal(m) = 0.0; pesystem(m) = 0.0; % initialize %
% Find n0 and the variance of the noise. %
string1 = ['Eb/No = ',num2str(ebn0db(m))];
disp(string1) % track execution 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); for n=1:nx
theta = angle(xx(n));
if (theta<0)
theta = theta+2*pi; end %
% Rotate x and y to the first quadrant and compute BER. %
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)); % reference d3 = real(yyy(n)); d4 = imag(yyy(n)); % system
pe1 = Qfunct(d1/sigma1)+Qfunct(d2/sigma1); % reference pe2 = Qfunct(d3/sigma)+Qfunct(d4/sigma); % system peideal(m) = peideal(m)+pe1; % SER of reference pesystem(m) = pesystem(m)+pe2; % SER of system end end
peideal = (1/2)*peideal./nx; % convert to BER pesystem = (1/2)*pesystem./nx; % conve
function [x, bits] = random_binary(nbits,nsamples)
% This function genrates a random binary waveform of length nbits % sampled 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.
function [ y ] = Qfunct( x ) y=(1/2)*erfc(x/sqrt(2)); end
%4QAM clear echo on
SNRindB1=0:2:15; SNRindB2=0:0.1:15; M=4;
k1=log2(M);
for i=1:length(SNRindB1),
smld_err_prb_4(i)=QAM_4(SNRindB1(i)); % 仿真误符号率
echo off; end; echo on ;
for i=1:length(SNRindB2),
SNR=exp(SNRindB2(i)*log(10)/10); %理论上误符号率
theo_err_prb_4(i)=4*Qfunct(sqrt(3*k1*SNR/(M-1))); echo off ; end; echo on ;
% 绘制误符号率对比图
semilogy(SNRindB1,smld_err_prb_4,'*'); hold on
semilogy(SNRindB2,theo_err_prb_4);
legend('4QAM仿真曲线','4QAM理论曲线',2);
function [p]=QAM_4(snr_in_dB) N=10000;
d=1; % 符号间最小距离 Eav=2*d^2; % 符号能量
snr=10^(snr_in_dB/10); % 每比特符号的信噪比 sgma=sqrt(Eav/(4*snr)); % 噪声方差 M=4;
% 产生数据 for i=1:N,
temp=rand;
dsource(i)=1+floor(M*temp); % 生成[1，4]的随机数 end;
%映射图案
mapping=[-d d; d d;
-d -d; d -d]; for i=1:N,
qam_sig(i,:)=mapping(dsource(i),:); end;
% 接收信号 for i=1:N,
[n(1) n(2)]=gngauss(sgma); r(i,:)=qam_sig(i,:)+n; end;
% 解调检测及误码率计算 numoferr=0; for i=1:N, for j=1:M,
metrics(j)=(r(i,1)-mapping(j,1))^2+(r(i,2)-mapping(j,2))^2; end;
[min_metric decis] = min(metrics); if (decis~=dsource(i)), numoferr=numoferr+1; end; end;
p=numoferr/(N);