QPSK通过Rayleigh信道多径衰落的Matlab仿真.doc

参照通信系统仿真原理与无线应用351页例14-1在这个例子里，我们对有3条固定路径的AWGN多径信道中的QPSK系统进行BER性能仿真，并与在理想的AWGN信道（没有多径）中同样系统地BER性能进行比较 书上有比较详细的数学推导，不抄了。这个例子似乎没有考虑多普勒频移。待我继续学习下一个例子，这个也没太看懂。 下面是该例子的源程序，P0、P1、P2分别是LOS路径和两条延迟瑞利分量的相对功率级。当p0=0且delay！=0时为瑞利频率选择性衰落，delay=0时为瑞利平坦衰落。 主程序scriptfile：% 两径瑞利衰落信道仿真% 设定默认参数NN=256; % 传输符号个数tb=0.5; % 一比特时间fs=10; % 每符号采样数ebn0db=1:2:15; % 设定Eb/N0% 建立QPSK信号x=random_binary(NN,fs)+i*random_binary(NN,fs); % x为QPSK信号% 输入功率和延迟p0=0; % 视距LOS分量p1=20; % 第一路径分量p2=1; % 第二路径分量delay=1; % 按照每符号采样数决定的延迟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); endendy1=x;for k=1:delay2 y2(1,k)=y1(1,k)*sqrt(p0);endfor 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);% 仿真结束% Use the semianalytic BER estimator . The following sets up the semi% analytic estimator . Find the maximun 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 calibrationhh=impz(b,a);ts=1/16;nbw=(fs/2)*sum(hh.2);% 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);figuresemilogy(ebn0db,peideal,b*-,ebn0db,pesystem,r+-)xlabel(Eb/N0 (db);ylabel(Probability of Error);grid on axis(0 14 10(-10) 1);% End of script file. 相关的一些调用程序（4个）：1 vxcorr.mfunction 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)=Ea(n+m)*conj(b(n)=Ea(n)*conj(b(n-m)a=a(:); % convert a to column vectorb=b(:); % convert b to column vectorM=length(a); % same as length(b)maxlag=M-1; % maximum value of laglags=-maxlag:maxlag;A=fft(a,2nextpow2(2*M-1); % fft of AB=fft(b,2nextpow2(2*M-1); % fft of Bc=ifft(A.*conj(B); % corsscorrelation% 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(nrnc) c=c.; lags=lags.;end% End of function file. 2 random_binary.mfunction x,bits=random_binary(nbits,nsamples)% This function generates 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); endend% End of function file. 3 qpsk_berest.m% File: psk_berest.mfunction peideal,pesystem=psk_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 filtered output;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 bandwidthfor 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);% % 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)