clc,clear,close;% Symbol sequence lengthL=100000;% Generate the original bit sequencesourceSeq=randnum(L);[pI,pQ,sourceSeqCode]=Map(sourceSeq,L);Eb = 1/3; errbit = zeros(1,26);errnum = zeros(1,26);SNR=-5:20;LS = length(SNR);for i=1:LS % Find the one-sided power spectral density of noise for a given signal-to-noise ratio N0=Eb/(10^(SNR(i)/10)); % variance var(i)=N0/2 ; [rI,rQ]=noise(var(i),pI,pQ); [result,I]=judgment(rI,rQ); [errbit(i),errnum(i)]=Count(result,I,sourceSeq,sourceSeqCode);enddraw(sourceSeqCode,rI,rQ);
复制代码
4.2 QPSK 与 8PSK 误码率对比主函数
1.流程图
图 4.2 qpsk 与 8psk 误码率对比主函数流程图
2.代码实现
clc,clear,close;% Symbol sequence lengthL=100000;% Generate the original bit sequencesourceSeq=randnum(L);[pI,pQ,sourceSeqCode]=Map(sourceSeq,L);Eb = 1/3; errbit = zeros(1,26);errnum = zeros(1,26);SNR=-5:20;LS = length(SNR);for i=1:LS % Find the one-sided power spectral density of noise for a given signal-to-noise ratio N0=Eb/(10^(SNR(i)/10)); % variance var(i)=N0/2 ; [rI,rQ]=noise(var(i),pI,pQ); [result,I]=judgment(rI,rQ); [errbit(i),errnum(i)]=Count(result,I,sourceSeq,sourceSeqCode);end% draw(sourceSeqCode,rI,rQ); % Calculate the theoretical bit error rate using the erfc functionTheory8PSKSER = erfc(sqrt(3*10.^(SNR/10)) * sin(pi/8)); % QPSK TheoryQPSKSER = erfc(sqrt(10.^(SNR/10))).*(1-0.25*erfc(sqrt(10.^(SNR/10))));;% Load the SER-SNR curve of QPSKload('qpsk_errnum');figure(2);epskerr = errnum/L;semilogy(SNR,epskerr,'r-o');hold on;semilogy(SNR,qpskerr,'b-o');hold on;semilogy(SNR,Theory8PSKSER,'k-');hold on;semilogy(SNR,TheoryQPSKSER,'k-');hold on;ylabel('SER');xlabel('SNR/dB')legend('8PSK仿真曲线','QPSK仿真曲线','理论曲线','Location', 'northeast' );grid on;axis([-5,20,10e-7,1]);
复制代码
五、子函数设计
5.1 随机比特序列的产生
代码实现:
function [SourceSeq]=randnum(L)% L is the length of the generated sequence code% Since a code is composed of 3 bits, it is generated here with 3*L.randnum=rand(3*L,1);% Initialize the original sequenceSourceSeq=zeros(3*L,1);% The randomly generated sequence is judged% if the random number is greater than 0.5, it is judged as 1, otherwise it is judged as 0.for i=1:3*L if(randnum(i)>=0.5) SourceSeq(i)=1; else SourceSeq(i)=0; endend
复制代码
5.2 格雷编码序列
代码实现
% Define 8psk mapping functionfunction [pI,pQ,SourceCode] = Map(SourceSeq,L)% pI - in-phase component% pQ - quadrature component% SourceCode - The size of the binary number of each digit of the sequence% initializationpI = zeros(L,1);pQ = zeros(L,1);% In order to facilitate subsequent expressions, sqrt(2)/2 is represented hereroot =sqrt(2)/2;% Constructing the mapping matrix according to the Gray code of 8PSKMappingMat = [[1,0];[root,root];[-root,root];[0,1];[root,-root];[0,-1];[-1,0];[-root,-root]]; SourceCode =zeros(L,1);% mapping processfor i=1:L % Since a source symbol is composed of three bits % the low and high bits are read in reverse order here, and expressed in decimal SourceCode(i)=SourceSeq(3*i-2)*4+SourceSeq(3*i-1)*2+SourceSeq(3*i)+1; % Find the corresponding code through the decimal representation and map it pI(i) = MappingMat(SourceCode(i),1); pQ(i) = MappingMat(SourceCode(i),2);endend % Octal Gray Code Conversionfunction[a1,a2]=Map_other(N)% a1 is a sequence of random bits in binary% a2 is a sequence of octal symbols% Random bit sequence of length 3La1=bit(3*N);% a2 is used to store the code element sequence of length La2=zeros(1,N);for i=1:3:3*N-2 a2((i+2)/3)=abs(a1(i)*7-abs(a1(i+1)*3-a1(i+2)));% Converts binary bit sequence Gray-encoded to octal sequenceend end
复制代码
5.3 映射函数
代码实现
% 8PSK coordinate mappingfunction [y3]=coordinate(x1,bit)% x1 is the encoded octal sequence, bit is the originally generated binary random sequenceN=length(x1);Es=bit*bit'/N;% The first line of y3 is used to store the abscissa% the second line is used to store the ordinatey3=zeros(2,N);% Coordinate mappingfor i=1:N y3(1,i)=sqrt(Es)*cos(pi/4*x1(i)+pi/8); y3(2,i)=sqrt(Es)*sin(pi/4*x1(i)+pi/8);endend
复制代码
5.4 噪声生成与叠加输出
代码实现
% Generate Gaussian random noise sub-function , var is the variancefunction [rI,rQ] = noise(var,pI,pQ)L = length(pI);nc=zeros(L,1);ns=zeros(L,1);for k=1:Lu=rand;z=sqrt(var*2*log(1/(1-u)));nc(k)=z*cos(2*pi*u);ns(k)=z*sin(2*pi*u);end% Output two mutually orthogonal Gaussian signalsrI = pI+nc;rQ = pQ+ns;end
复制代码
5.5 判决函数
代码实现
% Judgment Criterion: Minimum Euclidean Distance Criterionfunction [result,I]=judgment(rI,rQ)root=sqrt(2)/2;L = length(rI);I=zeros(L,1);mapl = [[1,0];[root,root];[-root,root];[0,1];[root,-root];[0,-1];[-1,0];[-root,-root]];for i=1:L index=0; minp=100; for j=1:8 % Traverse each coordinate % find the point with the smallest Euclidean distance from the point, and record it. if((mapl(j,1)-rI(i))^2+(mapl(j,2)-rQ(i))^2<minp) minp=(mapl(j,1)-rI(i))^2+(mapl(j,2)-rQ(i))^2; index=j; end end I(i)=index;end % According to the judgment point, make the corresponding assignmentresult = zeros(3*L,1);for i=1:L switch(I(i)) case 1 [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(0,0,0); case 2 [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(0,0,1); case 3 [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(0,1,0); case 4 [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(0,1,1); case 5 [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(1,0,0); case 6 [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(1,0,1); case 7 [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(1,1,0); case 8 [result(3*i-2),result(3*i-1),result(3*i)]=SetValue(1,1,1); endend end
function draw(I,rI,rQ)figure;root = sqrt(2)/2;%映射矩阵MappingMat = [[1,0];[root,root];[-root,root];[0,1];[root,-root];[0,-1];[-1,0];[-root,-root]];L=length(I)for i=1:L % Draw points with different colors according to different values switch(I(i)) case 1 plot(rI(i),rQ(i),'*','color','r');hold on; case 2 plot(rI(i),rQ(i),'*','color','g');hold on; case 3 plot(rI(i),rQ(i),'*','color','b');hold on; case 4 plot(rI(i),rQ(i),'*','color','c');hold on; case 5 plot(rI(i),rQ(i),'*','color','m');hold on; case 6 plot(rI(i),rQ(i),'*','color','y');hold on; case 7 plot(rI(i),rQ(i),'*','color','k');hold on; case 8 plot(rI(i),rQ(i),'*','color','[0.5,0.5,0.5]');hold on; endendx = -4:0.1:4;y = -2:0.1:2;% plot(x,zeros(length(x)),'k');hold on;% plot(zeros(length(y)),y,'k');hold on;axis equalaxis([-2,2,-2,2]);end
复制代码
5.7 误码率及误比特率计算函数
代码实现
function [errbit,errsymbol]=Count(result,I,sourceSeq,sourceSym)% result: the sequence after the decision% I: symbol after decision% sourceSeq: original bit sequence% sourceSym: original symbol sequenceerrbit=0;errsymbol=0;for i=1:length(sourceSeq) if(sourceSeq(i)~=result(i)) errbit=errbit+1; endendfor i=1:length(I) if(I(i)~=sourceSym(i)) errsymbol=errsymbol+1; endendend
复制代码
六、性能分析与实验结果
6.1 比较 8PSK 与 QPSK 的 Monte Carlo 仿真误符号率曲线、理论误符号率曲线差别
在 AWGN 信道下,未加信道纠错码的 8PSK 调制通信系统检测器的判决准则选为最小距离法(星座图上符号间的距离),格雷码映射,比较数据点为 100000 时 8PSK 与 QPSK 的 Monte Carlo 仿真误符号率曲线,理论误符号率曲线,比较差别(横坐标是 SNR=Eb/N0)。(一张图上呈现 4 条曲线)
图 6. 1 QPSK 与 8PSK 性能比较
通过 Monte Carlo 仿真误符号率曲线可以看出,整体仿真结果基本符合理论计算曲线,并且在不同的信噪比下,对应的 QPSK 的误码率明显小于 8PSK 的误码率,8PSK 的星座图如下:
评论