基于matlab的QPSK与BPSK信号性能比较仿真

┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊

目录

第一章概述 (1)

第二章QPSK通信系统原理与仿真 (1)

2.1 QPSK系统框图介绍 (1)

2.2QPSK信号的调制原理 (2)

2.2.1QPSK信号产生方法 (2)

2.2.2QPSK星座图 (2)

2.3QPSK解调原理及误码率分析 (3)

2.3.1QPSK解调方法 (3)

2.3.2QPSK系统误码率 (3)

2.4QPSK信号在AWGN信道下仿真 (4)

第三章BPSK通信系统原理与仿真 (4)

3.1BPSK信号的调制原理 (4)

3.2BPSK解调原理及误码率分析 (4)

第四章QPSK与BPSK性能比较 (5)

4.1QPSK与BPSK在多信道下比较仿真 (5)

4.1.1纵向比较分析 (5)

4.1.2横向比较分析 (7)

4.2仿真结果分析 (7)

4.2.1误码率分析 (7)

4.2.2频带利用率比较 (7)

附录 (8)

代码1 (8)

代码2 (8)

代码3 (10)

代码4 (12)

┊

┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊

第一章概述

QPSK是英文Quadrature Phase Shift Keying的缩略语简称，意为正交相移键控，是一种数字调制方式。它以其抗干扰性能强、误码性能好、频谱利用率高等优点，广泛应用于数字微波通信系统、数字卫星通信系统、宽带接人、移动通信及有线电视系统之中。

BPSK是英文Binary Phase Shift Keying的缩略语简称，意为二相相移键控，是利用偏离相位的复数波浪组合来表现信息键控移相方式的一种。它使用了基准的正弦波和相位反转的波浪，使一方为0，另一方为1，从而可以同时传送接受2值(1比特)的信息。

本文所研究的QPSK系统与二进制的BPSK系统相比，具有以下特点:

1.在传码率相同的情况下，四进制数字调制系统的信息速率是二进制系统的2倍。

2.在相同信息速率条件下，四进制数字调制系统的传码率是二进制系统的1/4倍，这一特

点使得四进制码元宽度是二进制码元宽度的2倍，码元宽度的加大，可增加每个码元的

能量，也可减小码间串扰的影响。

3.由于四进制码元速率比二进制的降低，所需信道带宽减小。

4.在接收系统输入信噪比相同的条件下，四进制数字调制系统的误码率要高于二进制系

统。

5.四进制数字调制系统较二进制系统复杂，常在信息速率要求较高的场合。

基于以上优点，在数字信号的调制方式中QPSK（Quadrature Phase Shift Keying）四相移键控是目前最常用的一种卫星数字信号调制方式,它具有较高的频谱利用率、较强的抗干扰性，在电路上实现也较为简单，因而被WCDMA和CDMA2000等第三代移动通信系统采用。

第二章QPSK通信系统原理与仿真

2.1 QPSK系统框图介绍

在图2.1的系统中，发送方，QPSK数据源采用随机生成，信源编码采用差分编码，编码后的信号经QPSK调制器，经由发送滤波器进入传输信道。

接收方，信号首先经过相位旋转，再经匹配滤波器解调，经阈值比较得到未解码的接收信号，差分译码后得到接收信号，与信源发送信号相比较，由此得到系统误码率，同时计算系统误码率的理论值，将系统值与理论值进行比较。

对于信道，这里选取的是加性高斯白噪声（Additive White Gaussian Noise）以及多径Rayleigh

┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊衰落信道（Multipath Rayleigh Fading Channel）。

图2.1 QPSK系统框图

在实验中，选用的是差分码。

差分码又称为相对码，在差分码中利用电平跳变来分别表示1或0，分为传号差分码和空号差分码。

传号差分码：当输入数据为“1”时，编码波型相对于前一码电平产生跳变；输入为“0”时，波型不产生跳变。

空号差分码：当输入数据为“0”时，编码波型相对于前一码电平产生跳变；输入为“1”时，波型不产生跳变。

2.2QPSK信号的调制原理

2.2.1QPSK信号产生方法

QPSK信号通常由图2.2所示的调相法产生，输入的二进制序列经过一个串/并转换器后分为两路二进制序列，这里假设两路序列为a，b，并通过平衡调制器分别对同相载波及正交载波进行二相调制，得到如图2.3中的虚线矢量。

QPSK信号每个码元含有2比特的信息，如图2.3所示现用ab代表这两个比特。两个比特有4种组合，即00、01、10和11。它们和相位之间的关系通常都按格雷码的规律安排，如表2-1所示。

图2.2 调相法产生QPSK信号图2.3 信号矢量

2.2.2QPSK星座图

在本次系统仿真中使用仿真代码产生QPSK信号星座图，具体代码如附录代码1。

┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊运行得到的信号星座图如图2.4，2.5：

表2-1 QPSK编码规则

表中

k

－一组间隔均匀的受调制相位

2.3QPSK解调原理及误码率分析

2.3.1QPSK解调方法

由于

QPSK信号可以看作两个正交2PSK信号的叠加，解调框图如图2.6，用相干解调方法，即用两路正交的相干载波，可以很容易的分离出这两路正交的2PSK信号。解调后的两路基带信号码元a和b，经过并串变换后，成为串行数据输出。

图2.4 A方式信号星座图图2.5 B方式信号星座图

2.3.2QPSK系统误码率

在QPSK体制中，由其矢量图（图2.7）可以看出，因噪声的影响使接收端解调时发生错误判决，是由于信号矢量的相位发生偏离造成的。例如，设发送矢量的相位为

45，它代表基带信号码元“11”，若因噪声的影响使接收矢量的相位变成

135，则将误判为“01”。当各个发送矢量以等概率出现时，合理的判决门限应该设在和相邻矢量等距离的位置。在图中对于矢量“11”来说，判决门限应该设在 0和

90。当发送“11”时，接收信号矢量的相位若超出这一范围（图图2.6 QPSK信号解调原理框图图2.7 QPSK噪声容限

┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊中阴影区），则将发生错判。设)

(θ

f为接收矢量（包括信号和噪声）相位的概率密度，则发生错误的概率为：

省略计算)

(θ

f和Pe的繁琐过程，直接给出计算结果：

上式计算出的是QPSK信号的误码率。若考虑其误比特率，正交的两路相干解调方法和2PSK 中采用的解调方法一样。所以其误比特率的计算公式也和2PSK的误码率公式一样。

2.4QPSK信号在AWGN信道下仿真

仿真代码见附录代码2。

运行结果如图2.8：

图2.8 AWGN环境下QPSK系统的仿真和理论结果

由运行结果可以看出，仿真值与理论值基本吻合，误码率值随信噪比增大而下降。

第三章BPSK通信系统原理与仿真

3.1BPSK信号的调制原理

图3.1 模拟调制方法

图3.2 键控调制方法

BSPK信号通常有2种调制方式，分别如图3.1、3.2所示。

在2PSK中，通常用初始相位0和p分别表示二进制“1”和“0”。因此，2PSK信号的时域表达式为

式中，

n

?表示第n个符号的绝对相位：

因此，上式可以改写为

3.2BPSK解调原理及误码率分析

BPSK信号的解调器原理方框图3.3和波形图3.4，这里给出的解调框图采用相干解调的接收方法。

图3.3 BPSK信号相干解调框图

波形图中，假设相干载波的基准相位与2PSK信号的调制载波的基准相位一致（通常默认为0相位）。但是，由于在2PSK信号的载波恢复过程中存在着的相位模糊，即恢复的本地载波与

┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊所需的相干载波可能同相，也可能反相，这种相位关系的不确定性将会造成解调出的数字基带信号与发送的数字基带信号正好相反，即“1”变为“0”，“0”变为“1”，判决器输出数字信

图3.4 相干解调波形示意图

号全部出错。这种现象称为2PSK 方式的“倒π”现象或“反相工作”。这也是2PSK方式在实际中很少采用的主要原因。另外，在随机信号码元序列中，信号波形有可能出现长时间连续的正弦波形，致使在接收端无法辨认信号码元的起止时刻。

第四章QPSK与BPSK性能比较

4.1QPSK与BPSK在多信道下比较仿真

QPSK仿真程序见代码3。BPSK仿真程序见代码4。

下面对QPSK与BPSK在多信道下的BER-SNR曲线进行横向及纵向的比较，如表4-1所示。

4.1.1纵向比较分析

图中依次给出LOS信道、莱斯平坦衰落信道、莱斯频率选择性衰落、瑞利平坦衰落信道、瑞利频率选择性衰落。比较频率选择性衰落和其相应的平坦衰落，可发现频率选择性衰落比相应的平坦衰落的性能差。比较莱斯信道和瑞利信道，可发现不存在视距分量的瑞利信道比存在视距分量的任何莱斯信道的性能差。

表4-1 QPSK与BPSK在多信道下性能比较

┊

┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊

4.1.2横向比较分析

相同信道下，BPSK调制的系统误码率小于QPSK调制。本次仿真过程中，莱斯平坦衰落信道下两者误码率差别显著。

4.2仿真结果分析

4.2.1误码率分析

由前面介绍的误码率内容可以看出，QPSK判决门限为

90，BPSK的判决门限为

180。因此相同系统情况下的误码率BPSK优于QPSK。

4.2.2频带利用率比较

在传码率相同的情况下，四进制数字调制系统的信息速率是二进制系统的2倍。频带利用率公式

在相同信号速率的情况下，QPSK和BPSK系统的带宽是相同的，但是由于QPSK每个信号都是四进制的，QPSK每个信号包含2bit信息，所以比特率就是BPSK的两倍，因而其频带利用率即为BPSK的两倍。BPSK系统理论的频带利用率最大为1，但是在实际的实现中不能达到1，而在QPSK系统中，频带利用率可以超过1。

┊┊┊

┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊

附录

代码1

function qpskconstellationA(M)

M=4;

x=[0:M-1];

scatterplot(pskmod(x,M)); % A方式QPSK信号

grid on;

return

function qpskconstellationB(M)

M=4;

x=[0:M-1];

scatterplot(pskmod(x,M,pi/4)); % B方式QPSK信号

grid on;

return

代码2

function [x,bits]=random_binary(nbits,nsamples)

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

return

function [c,lags]=vxcorr(a,b)

%计算a b 的互相关系数，返回C是列向量，长度为ab的长度和减1

a=a(:);

b=b(:);

M=length(a);

maxlag=M-1;

lags=[-maxlag:maxlag]';

A=fft(a,2^nextpow2(2*M-1));%快速傅立叶变换

B=fft(b,2^nextpow2(2*M-1));

c=ifft(A.*conj(B));%快速傅立叶反变换

c=[c(end-maxlag+1:end,1);c(1:maxlag+1,1)];

[nr nc]=size(a);

if(nr>nc)

c=c';

lags=lags.';

end

return

function BER_MC=MCQPSKrun(N,Eb,N0,ChanAtt,TimingBias,TimingJitter,PhaseBias,PhaseJitter) fs=1e+6; % 采样速率

SymRate=1e+5; % 信号速率

Ts=1/fs; % 采样周期

TSym=1/SymRate; % 信号周期

SymToSend=N; % 发送信号数

ChanBW=4.99e+5; % 带宽

MeanCarrierPhaseError=PhaseBias; % 载波相位均值StdCarrierPhaseError=PhaseJitter; % 相位误差标准差MeanSymbolSyncError=TimingBias; % 符号同步误差均值

┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊

StdSymbolSyncError=TimingJitter; % 符号同步误差标准差

ChanGain=10^(-ChanAtt/20); % 信道增益

TxBitClock=Ts/2; % 发送机时钟

RxBitClock=Ts/2; % 接收机时钟

%接收机输入端噪声标准差和信号幅度

RxNoiseStd=sqrt((10^((N0-30)/10))*(fs/2));

TxSigAmp=sqrt(10^((Eb-30)/10)*SymRate);

%分配缓存

SampPerSym=fs/SymRate;

probe1=zeros((SymToSend+1)*SampPerSym,1);

probe1counter=1;

probe2=zeros((SymToSend+1)*SampPerSym,1);

probe2counter=1;

%计已传输信号的个数

TxSymSent=1;

RxSymDemod=0;

%发送和接收数据缓冲区

[unused,SourceBitsI]=random_binary(SymToSend,1);

[unused,SourceBitsQ]=random_binary(SymToSend,1);

%差分编码

TxBitsI=SourceBitsI*0;

TxBitsQ=SourceBitsQ*0;

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),So urceBitsQ(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

%产生复信号

TxBits=((TxBitsI*2)-1)+(sqrt(-1)*((TxBitsQ*2)-1));

RxIntegrator=0; %初始化接收机积分器

TxBitClock=2*TSym; %初始化发送机

%设计信道滤波器，产生滤波器参数序列

[b,a]=butter(2,ChanBW/(fs/2));

b=[1];

a=[1];

[junk,FilterState]=filter(b,a,0);

% 开始仿真循环

while TxSymSent

TxBitClock=TxBitClock+Ts;

if TxBitClock>TSym

TxSymSent=TxSymSent+1;

TxBitClock=mod(TxBitClock,TSym);

TxOutput=TxBits(TxSymSent)*TxSigAmp;

end

%信号经过信道滤波器

[Rx,FilterState]=filter(b,a,TxOutput,FilterState);

%加高斯白噪声

Rx=(ChanGain*Rx)+(RxNoiseStd*(randn(1,1)+sqrt(-1)*randn(1,1)));

%基于接收机载波同步误差的相位旋转

PhaseRotation=exp(sqrt(-1)*2*pi*(MeanCarrierPhaseError+(randn(1,1)*StdCarrierPhaseError))/360);

Rx=Rx*PhaseRotation;

probe1(probe1counter)=Rx;

probe1counter= probe1counter+1;

%更新接收机积分清除器

RxIntegrator=RxIntegrator+Rx;

probe2(probe2counter)=RxIntegrator;

probe2counter= probe2counter+1;

%更新接收机时钟，判断是不是适合采样

RxBitClock=RxBitClock+Ts;

RxTSym=TSym*(1+MeanSymbolSyncError+(StdSymbolSyncError*randn(1,1)));

if RxBitClock>RxTSym %解调信号

RxSymDemod=RxSymDemod+1;

┊┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊

RxBitsI(RxSymDemod)=round(sign(real(RxIntegrator))+1)/2;

RxBitsQ(RxSymDemod)=round(sign(imag(RxIntegrator))+1)/2;

RxBitClock=RxBitClock-TSym;

RxIntegrator=0;

end

end

%差分解码

SinkBitsI=SourceBitsI*0;

SinkBitsQ=SourceBitsQ*0;

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

%在输入和输出100字节中寻找最佳时延

[C,Lags]=vxcorr(SourceBitsI(10:110),SinkBitsI(10:110));

[MaxC,LocMaxC]=max(C);

BestLag=Lags(LocMaxC);

% 调整时延

if BestLag>0

SourceBitsI=SourceBitsI(BestLag+1:length(SourceBitsI));

SourceBitsQ=SourceBitsQ(BestLag+1:length(SourceBitsQ));

elseif BestLag<0

SinkBitsI=SinkBitsI(-BestLag+1:length(SinkBitsI));

SinkBitsQ=SinkBitsQ(-BestLag+1:length(SinkBitsQ));

end

%将序列调整成相同长度

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);

Errors=sum(SourceBitsI ~= SinkBitsI)+sum(SourceBitsQ ~= SinkBitsQ);

BER_MC=Errors/(2*length(SourceBitsI));

return

function MCQPSKBER

Eb=22:0.5:26;

N0=-50;

ChannelAttenuation=70;

EbN0dB=(Eb-ChannelAttenuation)-N0;

EbN0=10.^(EbN0dB./10);

BER_T=0.5*erfc(sqrt(EbN0));

N=round(100./BER_T);

BER_MC=zeros(size(Eb));

for k=1:length(Eb)

BER_MC(k)=MCQPSKrun(N(k),Eb(k),N0,ChannelAttenuation,0,0,0,0);

disp(['simulation',num2str(k*100/length(Eb)),'%complate']);

end

figure(1);

semilogy(EbN0dB,BER_MC,'o',EbN0dB,2*BER_T,'-');

xlabel('Eb/N0(dB)');

ylabel('bit error rate');

legend('MC BER Estimate','Theoretical BER');grid;

return

代码3

function [peideal,pesystem]=qpsk_berest(xx,yy,ebn0db,eb,tb,nbw)

%xx为输入序列，yy为输出序列，ebn0db为信噪比，eb为单个信号能力，tb信号周期，nbw噪声带宽

[n1 n2]=size(xx);

nx=n1*n2;

[n3 n4]=size(yy);

┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊[n5 n6]=size(ebn0db);

neb=n5*n6;

%接收机带宽设定为rs/2

nbwideal=1/(2*tb*2);

for m=1:neb

peideal(m)=0.0; % 初始化

pesystem(m)=0.0; % 初始化

%计算n0和噪声方差

string1=['Eb/N0=',num2str(ebn0db(m))];

disp(string1)

ebn0(m)=10^(ebn0db(m)/10);

n0=eb/ebn0(m); % 噪声功率

sigma=sqrt(n0*nbw*2); % 方差

sigma1=sqrt(n0*nbwideal*2); % 理想方差

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

%接收信号的相位旋转

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));

d3=real(yyy(n));

d4=imag(yyy(n));

pe1=q(d1/sigma1)+q(d2/sigma1);

pe2=q(d3/sigma)+q(d4/sigma);

peideal(m)=peideal(m)+pe1;

pesystem(m)=pesystem(m)+pe2;

end

end

peideal=(1/2)*peideal./nx;

pesystem=(1/2)*pesystem./nx;

return

function QPSKThreeray(p0,p1,p2,delay)

NN=256; %发送的信号数

tb=0.5; %信号周期

fs=16; %每信号的采样数

ebn0db=[1:2:14]; %信噪比矢量

%产生QPSK信号

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

delay0=0;

delay1=0;

delay2=delay;

%设置瑞利信道

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

%匹配滤波器

b=-ones(1,fs);

┊┊┊┊┊┊┊┊┊┊┊装┊┊┊┊┊订┊┊┊┊┊线┊┊┊┊┊┊┊┊┊┊┊┊┊a=1;

y=filter(b,a,y2);

%应用半解析方法，进行误码率估计。

[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);

hh=impz(b,a);

ts=1/fs;

nbw=(fs/2)*sum(hh.^2);

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');

ylabel('Probability of Error');

grid

axis([0 14 10^(-10) 1])

return

代码4

function bpskThreeray(p0,p1,p2,delay)

NN=256;

tb=0.5;

fs=16;

ebn0db=[1:2:14];

x=random_binary(NN,fs);

delay0=0;

delay1=0;

delay2=delay;

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

b=-ones(1,fs);

b=b/fs;

a=1;

y=filter(b,a,y2);

[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);

hh=impz(b,a);

ts=1/fs;

nbw=(fs/2)*sum(hh.^2);

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');

ylabel('Probability of Error');

grid

axis([0 14 10^(-10) 1])

return