MATLAB_QPSK调制与解调-(4820)

实验名称： QPSK仿真系

统一、实验目的：

1、学会 QPSK 调制与解调系统的构成

2、学会 QPSK 调制与解调系统的各模块的构建

3、学会误码率与误符号率的统计方法以及Matlab 算法

二、实验原理：

1、QPSK：四进制绝对相移键控，也称为多进制数字相位调制，利用载波的四种

不同相位状态来表征数字信息的调制方式。

2、QPSK 的调制方法有正交调制方式（双路二相调制合成法或直接调相法）、相位选择法、插入脉冲法。

调制与解调系统的构成：

AWGN

产生器

n

c

s r

I I

b?b

0,1序列QPSK

r QPSK

s判决器产生器映射Q反映射

Q

n

s

AWGN

产生器

比较器

误比特统计误符号统计

3、各模块的实现方法：

（1）、信源的产生：使用randint(m,n,2)函数产生一个 m 行 n 列的随机二进制数列

（2）、QPSK 符号映射：将产生的 0， 1 比特流按照 QPSK 调制方式进行映射，本实验采用π/4 QPSK 的调制方式，图为：

1000

0111

（3）、AWGN 信号产生： AWGN 产生器就是产生满足均值为 0，方差为 1 的高斯白噪声。实验中使用 randn(m,n)函数产生一个 m 行 n 列的高斯噪声序列。

（4）、信号幅度控制：根据 AWGN 信道模型，接收信号可以分别表示为

r Q s Q n Q

α 就是当噪声功率归一化为r I s I n I

1（0均值，方差为 1）时，根据信噪比关系而计算

出来的信号平均幅度

SNR 10log v s2v

s2

SNR v2sqrt v n10 ^10 n

v s

（ 5）、QPSK 反映射及判决：对接收到的信号在 4 种可能的四种信号向量 [(1,0), (0,1), (-1,0), (0,-1)] 上投影 (即进行点积 )。投影最大的值所对应的信号向量就是所发送信号的符号值，然后恢复出比特流

（6）、误码率及误符号率统计：

误码率：将检测出来的比特流和发送的原始比特流进行比较，统计出出现错误的比特数

误符号率：将检测出来的比特流变成两组，构成符号，和发送端符号映射后的符

号流进行比较，只要符号中任错一 bit，就算该符号出错。统计出现错误的符号

数

三、实验内容：

1、调制与解调

clear all

close all

%调制

bit_in = randint(1e3, 1, [0 1]);

bit_I= bit_in(1:2:1e3);

bit_Q= bit_in(2:2:1e3);

data_I = -2*bit_I+1;

data_Q = -2*bit_Q+1;

data_I1=repmat(data_I',20,1);

data_Q1=repmat(data_Q',20,1);

for i=1:1e4

data_I2(i)=data_I1(i);

data_Q2(i)=data_Q1(i);

end;

f=0:0.1:1;

xrc=0.5+0.5*cos(pi*f);

data_I2_rc=conv(data_I2,xrc)/5.5;

data_Q2_rc=conv(data_Q2,xrc)/5.5;

f1=1;

t1=0:0.1:1e3+0.9;

n0=rand(size(t1));

I_rc=data_I2_rc.*cos(2*pi*f1*t1);

Q_rc=data_Q2_rc.*sin(2*pi*f1*t1); QPSK_rc=(sqrt(1/2).*I_rc+sqrt(1/2).*Q_rc); QPSK_rc_n0=QPSK_rc+n0;

%解调

I_demo=QPSK_rc_n0.*cos(2*pi*f1*t1);

Q_demo=QPSK_rc_n0.*sin(2*pi*f1*t1);

%低通滤波

I_recover=conv(I_demo,xrc);

Q_recover=conv(Q_demo,xrc);

I=I_recover(11:10010);

Q=Q_recover(11:10010);

t2=0:0.05:1e3-0.05;

t3=0:0.1:1e3-0.1;

%抽样判决

data_recover=[];

for i=1:20:10000

data_recover=[data_recover I(i:1:i+19) Q(i:1:i+19)];

end;

bit_recover=[];

for i=1:20:20000

if sum(data_recover(i:i+19))>0

data_recover_a(i:i+19)=1;

bit_recover=[bit_recover 1];

else

data_recover_a(i:i+19)=-1;

bit_recover=[bit_recover -1];

end

end

error=0;

dd= -2*bit_in+1;

ddd=[dd'];

ddd1=repmat(ddd,20,1);

for i=1:2e4

ddd2(i)=ddd1(i);

end

for i=1:1e3

if bit_recover(i)~=ddd(i)

error=error+1;

end

end

p=error/1000;

figure(1)

subplot(2,1,1);plot(t2,ddd2);axis([0 100 -2 2]);title(' 原序列 ');

subplot(2,1,2);plot(t2,data_recover_a);axis([0 100 -2 2]);title(' 解调后序列 ');

2、误码率仿真

%QPSK 误码率分析

SNRindB1=0:2:10;

SNRindB2=0:0.1:10;

for i=1:length(SNRindB1)

[pb,ps]=cm_sm32(SNRindB1(i));

smld_bit_err_prb(i)=pb;

smld_symbol_err_prb(i)=ps;

end;

for i=1:length(SNRindB2)

SNR=exp(SNRindB2(i)*log(10)/10);

theo_err_prb(i)=Qfunct(sqrt(2*SNR)); end;

title('QPSK 误码率分析 ');

semilogy(SNRindB1,smld_bit_err_prb,'*');

axis([0 10 10e-8 1]);

hold on;

%semilogy(SNRindB1,smld_symbol_err_prb,'o'); semilogy(SNRindB2,theo_err_prb); legend('仿真比特误码率 ','理论比特误码率 '); hold off;

function[y]=Qfunct(x)

y=(1/2)*erfc(x/sqrt(2));

function[pb,ps]=cm_sm32(SNRindB)

N=10000;

E=1;

SNR=10^(SNRindB/10);

sgma=sqrt(E/SNR)/2;

s00=[1 0];

s01=[0 1];

s11=[-1 0];

s10=[0 -1];

for i=1:N

temp=rand;

if (temp<0.25)

dsource1(i)=0;

dsource2(i)=0;

elseif (temp<0.5)

dsource1(i)=0;

dsource2(i)=1;

elseif (temp<0.75)

dsource1(i)=1;

dsource2(i)=0;

else

dsource1(i)=1;

dsource2(i)=1;

end;

end;

numofsymbolerror=0;

numofbiterror=0;

for i=1:N

n=sgma*randn(size(s00));

if((dsource1(i)==0)&(dsource2(i)==0))

r=s00+n;

elseif((dsource1(i)==0)&(dsource2(i)==1)) r=s01+n;

elseif((dsource1(i)==1)&(dsource2(i)==0)) r=s10+n;

else

r=s11+n;

end;

c00=dot(r,s00);

c01=dot(r,s01);

c10=dot(r,s10);

c11=dot(r,s11);

c_max=max([c00 c01 c10 c11]);

if (c00==c_max)

decis1=0;decis2=0;

elseif(c01==c_max)

decis1=0;decis2=1;

elseif(c10==c_max)

decis1=1;decis2=0;

else

decis1=1;decis2=1;

end;

symbolerror=0;

if(decis1~=dsource1(i))

numofbiterror=numofbiterror+1;

symbolerror=1;

end;

if(decis2~=dsource2(i))

numofbiterror=numofbiterror+1;

symbolerror=1;

end;

if(symbolerror==1)

numofsymbolerror=numofsymbolerror+1;

end;

end;

ps=numofsymbolerror/N;

pb=numofbiterror/(2*N);

3、QPSK 在 AWGN 信道下的仿真close all

clc

clear all

SNR_DB=[0:1:12];

sum=1000000;

data= randsrc(sum,2,[0 1]);

[a1,b1]=find(data(:,1)==0&data(:,2)==0);

message(a1)=-1-j;

[a2,b2]=find(data(:,1)==0&data(:,2)==1);

message(a2)=-1+j;

[a3,b3]=find(data(:,1)==1&data(:,2)==0);

message(a3)=1-j;

[a4,b4]=find(data(:,1)==1&data(:,2)==1);

message(a4)=1+j;

scatterplot(message)

title('B 点信号的星座图')

A=1;

Tb=1;

Eb=A*A*Tb;

P_signal=Eb/Tb;

NO=Eb./(10.^(SNR_DB/10));

P_noise=P_signal*NO;

sigma=sqrt(P_noise);

for Eb_NO_id=1:length(sigma)

noise1=sigma(Eb_NO_id)*randn(1,sum);

noise2=sigma(Eb_NO_id)*randn(1,sum);

receive=message+noise1+noise2*j;

resum=0;

total=0;

m1=find(angle(receive)<=pi/2&angle(receive)>0);

remessage(1,m1)=1+j;

redata(m1,1)=1;

redata(m1,2)=1;

m2= find( angle(receive)>pi/2&angle(receive)<=pi);

remessage(1,m2)=-1+j;

redata(m2,1)=0;

redata(m2,2)=1;

m3=find( angle(receive)>-pi&angle(receive)<=-pi/2);

remessage(1,m3)=-1-j;

redata(m3,1)=0;

redata(m3,2)=0;

m4=find( angle(receive)>-pi/2&angle(receive)<=0);

remessage(1,m4)=1-j;

redata(m4,1)=1;

redata(m4,2)=0;

[resum,ratio1]=symerr(data,redata);

pbit(Eb_NO_id)=resum/(sum*2);

[total,ratio2]=symerr(message,remessage);

pe(Eb_NO_id)=total/sum;

end

scatterplot(receive)

title('C 点信号的星座图')

Pe=1-(1-1/2*erfc(sqrt(10.^(SNR_DB/10)/2))).^2;

Pbit=1/2*erfc(sqrt(10.^(SNR_DB/10)/2));

figure(3)

semilogy(SNR_DB,pe,':s',SNR_DB,Pe,'-*',SNR_DB,pbit,'-o',SNR_DB,Pbit,':+')

legend('QPSK 仿真误码率 ','QPSK 理论误码率 ','QPSK 仿真误比特率','QPSK 理论误比特率',1) xlabel(' 信噪比 /dB')

ylabel(' 概率 P')

gird on

结果：

4、QPSK 在 AWGN 信道下的性能检测

clear all ;

%对随机序列产生的方法初始化

s = RandStream( 'mt19937ar' , 'Seed' , 5489);

RandStream.setDefaultStream(s);

s = RandStream( 'mcg16807' , 'Seed' ,0);

RandStream.setDefaultStream(s);

N = 10^5;%码元数

Eb_N0 = -4:20;%Eb/N0比特信噪比

s_data = zeros(1,N);

ErrorCount=zeros(1,length(Eb_N0));%预先分配内存，用于记录错误的码元数for j = 1:length(Eb_N0)

p = (2*(rand(1,N)>0.5)-1) + 1i*(2*(rand(1,N)>0.5)-1);

s = (1/sqrt(2))*p;%归一化

n = 10^(-Eb_N0(j)/40)*1/sqrt(2)*(randn(1,N) + 1i*randn(1,N));%加性高斯白噪声

r = s + n;

%解调

r_re = real(r);%实部

r_im = imag(r);%虚部

s_data(r_re < 0 & r_im < 0) = -1 + -1*1i;

s_data(r_re >= 0 & r_im > 0) = 1 + 1*1i;

s_data(r_re < 0 & r_im >= 0) = -1 + 1*1i;

s_data(r_re >= 0 & r_im < 0) = 1 - 1*1i;

ErrorCount(j) = size(find(p- s_data),2);

%错误的码元数

end

sim_QPSK = ErrorCount/N;

theory_QPSK = erfc(sqrt(0.5*(10.^(Eb_N0/20)))) -

(1/4)*(erfc(sqrt(0.5*(10.^(Eb_N0/20))))).^2;

close all

figure

semilogy(Eb_N0,theory_QPSK, 'b.-' );

hold on

semilogy(Eb_N0,sim_QPSK,'mx-' );

axis([-4 20 10^-3 1])

grid on