qam256.m

本文件是有关于QPSK仿真的matlab源码。

%	nnnn=1024;
%	ff=-nnnn/2:nnnn/2-1;
%	ff=ff/nnnn*2*fs/2;
%        plot(ff,data_fft);
%	title(' I Data Spectrum After Demodulation');
%	grid;
%	%pause;

%     low pass filtering

noise_i=filter(hrc_r,1,noise_i)*2;
noise_q=filter(hrc_r,1,noise_q)*2;

hrc_ra=zeros(size(hrc_r));
hrc_ra(81-ny_n:81+ny_n)=hrc_r(81-ny_n:81+ny_n);
%hrc_r=hrc_ra;
%????
if ny_full==1
	hrc_ra=hrc_r;
end;

rf_i_df=filter(hrc_ra,1,rf_i_d)*2;
rf_q_df=filter(hrc_ra,1,rf_q_d)*2;


%***************************
%	Add Channel Noise  *
%***************************

snr_analog=10*log10(abs(mean(rf_i_df.^2)/mean(noise_i.^2)));

snr_analog

rf_i_df=rf_i_df+noise_i;
rf_q_df=rf_q_df+noise_q;

%for itt=1:15

%	plot(rf_i_df(1+40*(itt-1):1+40*itt),'r*')
%	hold on;
%	plot(rf_i_df(1+40*(itt-1):1+40*itt))
%	hold off;
%	title(' I channel signal');
%	grid;
%	%pause;
%end




%       plot spectrum;
%	data_fft=20*log10(abs(fftshift(fft(kaiser(nnnn,11.)'.*...
%		 rf_i_df(offset:offset+1024-1)))));
%	nnnn=1024;
%	ff=-nnnn/2:nnnn/2-1;
%	ff=ff/nnnn*2*fs/2;
%       plot(ff,data_fft);
%	title(' I Data Spectrum After Demodulation Filter');
%	grid;
%	%pause;

clear rf_i_d rf_q_d rf_data if_data;
 
idat_ff=reshape(rf_i_df,4,length(rf_i_df)/4);
qdat_ff=reshape(rf_q_df,4,length(rf_i_df)/4);

%*********************************************
%	Demodulated I/Q Data: idat_f, qdat_f *
%*********************************************

% Sampling:
%	idat_f, qdat_f is the sampled output before the 
%	equalization;

idat_f=idat_ff(1,:);
qdat_f=qdat_ff(1,:);


clear idat_ff qdat_ff;

iq_f=idat_f+j*qdat_f;
plot(iq_f(200:800),'*');
title('Demod I/Q Constellation before equalization');
%axis([0 4 0 4]);
grid;
%pause

% Setup Equalizer

setup_a=0;
a=[0 0 0 0 0 ];
a=[a a a a ];
b=a;
a(13)=1;		% y(11) center tap =1;
y=a+1i*b;
y=y';			% y is the FFE coefficient;
e_length=20;		% equalizer length = 20 taps;

grp_delay=40+1;
grp_delay_w=n_test; %floor(e_length/2);	%for FBE
iq_w_s=iq_f*0;
off_set_n=30;

max_iie=length(iq_f)-grp_delay-e_length;
for iie=1:end_iie; %max_iie  % start the equalization loop

if round(iie/50)*50-iie==0
	plot(real(y),'*r');
	hold on
	plot(imag(y),'*g');
	title('Tap Value (imag(y), green)');
   grid
   axis([0 20 -.2 .2]);
   hold off
	%%pause
end
%????

%x=iq_f(iie+grp_delay+ch_delay:iie+grp_delay+e_length-1+ch_delay);	% delay line buffer

	% with training:

if training==1

	x(1:n_test)=idat(iie:iie+n_test-1)...
	       +1i*qdat(iie:iie+n_test-1);	% delay line buffer
	% no training:

else
	x(1:n_test)=iq_w_s(iie-n_test+off_set_n:iie-1+off_set_n);	% delay line buffer
end %of if

x(n_test+1:e_length)=iq_f(iie+grp_delay+ch_delay+n_test:iie+grp_delay+ch_delay+e_length-1);	% delay line buffer

x_i=idat(iie+grp_delay_w); % reference data;
x_q=qdat(iie+grp_delay_w); % reference data;

x_ii=idat(iie:iie+e_length-1); % reference data;
x_qq=qdat(iie:iie+e_length-1); % reference data;

%plot(x_ii,'r*');
%hold on
%plot(real(x),'g+');
%hold off
%grid
%title(' i data and reference in red');
%%pause

%plot(x_q,'r*');
%hold on
%plot(imag(x),'g*');
%hold off
%grid
%title(' q data and reference in red');
%%pause


w=x*y; 			% complex equalization;
w_i=real(w);
w_q=imag(w);

% 4 bit slicer:

% 	Thresholding 
if w_i > 8
	w_i=8;
        elseif w_i < -8
	w_i=-8;
end

if w_q > 8
	w_q=8;
        elseif w_q < -8
	w_q=-8;
end

% slicing

if abs(w_i-7.5)<=0.5
	w_i_s=7.5;
	elseif abs(w_i-6.5)<=0.5
	w_i_s=6.5;
	elseif abs(w_i-5.5)<=0.5
	w_i_s=5.5;
	elseif abs(w_i-4.5)<=0.5
	w_i_s=4.5;
	elseif abs(w_i+4.5)<=0.5
	w_i_s=-4.5;
	elseif abs(w_i+5.5)<=0.5
	w_i_s=-5.5;
	elseif abs(w_i+6.5)<=0.5
	w_i_s=-6.5;
	elseif abs(w_i+7.5)<=0.5
      w_i_s=-7.5;
  	elseif abs(w_i-3.5)<=0.5
	w_i_s=3.5;
	elseif abs(w_i-2.5)<=0.5
	w_i_s=2.5;
	elseif abs(w_i-1.5)<=0.5
	w_i_s=1.5;
	elseif abs(w_i-.5)<=0.5
   w_i_s=.5;
  	elseif abs(w_i+.5)<=0.5
	w_i_s=-.5;
	elseif abs(w_i+1.5)<=0.5
	w_i_s=-1.5;
	elseif abs(w_i+2.5)<=0.5
	w_i_s=-2.5;
	elseif abs(w_i+3.5)<=0.5
	w_i_s=-3.5;

end


if abs(w_q-7.5)<=0.5
	w_q_s=7.5;
	elseif abs(w_q-6.5)<=0.5
	w_q_s=6.5;
	elseif abs(w_q-5.5)<=0.5
	w_q_s=5.5;
	elseif abs(w_q-4.5)<=0.5
	w_q_s=4.5;
	elseif abs(w_q+4.5)<=0.5
	w_q_s=-4.5;
	elseif abs(w_q+5.5)<=0.5
	w_q_s=-5.5;
	elseif abs(w_q+6.5)<=0.5
	w_q_s=-6.5;
	elseif abs(w_q+7.5)<=0.5
      w_q_s=-7.5;
  	elseif abs(w_q-3.5)<=0.5
	w_q_s=3.5;
	elseif abs(w_q-2.5)<=0.5
	w_q_s=2.5;
	elseif abs(w_q-1.5)<=0.5
	w_q_s=1.5;
	elseif abs(w_q-.5)<=0.5
   w_q_s=.5;
  	elseif abs(w_q+.5)<=0.5
	w_q_s=-.5;
	elseif abs(w_q+1.5)<=0.5
	w_q_s=-1.5;
	elseif abs(w_q+2.5)<=0.5
	w_q_s=-2.5;
	elseif abs(w_q+3.5)<=0.5
	w_q_s=-3.5;

end



w_s=w_i_s+1i*w_q_s;

%plot(w,'r*');
%hold on
%plot(w_s,'*');
%hold off;
%title('w red and w_s');
%axis([-8 8 -8 8]);
%%pause
%	with training

if training==1;
	error_out=w-x_i-1i*x_q;    	% with training
	error_out_t=w-x_i-1i*x_q; 
	else
%	no training;
	error_out=w-w_s;		% no training;
        error_out_t=w-x_i-1i*x_q; 
	bit_error_b=abs(real(w_s-x_i-1i*x_q))+ 1i*abs(imag(w_s-x_i-1i*x_q)); 	% output bit error;
end;

%plot(x_i, '*r');
%hold on;
%plot(real(w_s), '+');
%hold off;
%grid
%%pause;

error_in(iie)=abs(x(n_test+1)-x_i-1i*x_q); %  error before equalization
error_plot(iie)=abs(error_out);		% observed error after equalization
error_out_tt(iie)=abs(error_out_t);	% true error after equalization.
bit_error(iie)=bit_error_b;		% bit error

%tt_iie=200;
%if round(iie/tt_iie)*tt_iie-iie==0
%   [ %y(test_tap1); ...
%      %x(test_tap1); ...
%      abs((y(test_tap1:test_tap2).*x(test_tap1:test_tap2)')); ...
%   %sum(abs(y(test_tap1:test_tap2))); ...
%   abs( error_out)]
%   iie
%   %pause
%end

%????


%	y(test_tap1:test_tap2)=y(test_tap1:test_tap2)-mu*error_out*(x(test_tap1:test_tap2)');
	y(test_tap1:test_tap2)=y(test_tap1:test_tap2)-mu*error_out*sign(x(test_tap1:test_tap2)');

iq_data_e(iie)=w;		% output equalization results
iq_w_s(iie+off_set_n)=w_s;
iq_data_in(iie)=x_i+1i*x_q;

%********************************************************
%		Adaptive Control			*
%********************************************************

test_error=mean( error_plot(max([iie-40,1]):iie));	% mean over 21 points;

 if converge_i==0 | converge_i==2;		% no more adjustment after converge and  error <0.2;

if converge_i ~=2 &iie >= i_module 	& test_error <.3 %  
	test_tap1=13-5;
	test_tap2=13+5;
  mu=.005
   converge_i=2;
elseif iie >= i_module 	& test_error < 0.2 %  
	test_tap1=1;
	test_tap2=19;
	mu=0.001
	converge_i=iie;
end
end


end 		% of iie;

plot(error_in,'r')
grid;
hold on;
plot(error_out_tt,'g');
plot(error_plot);
title(' error in (red) and error out observed (y) and true error(g)');

hold off;
%beep;
pause;

plot(real(bit_error))
grid
title('Symbol Error Plot, I channel');
pause;

%plot(error_plot(1:300));
%grid;
%title('error output');
%axis([1 iie-1 0 .1]);
%%pause

plot(real(y),'r*');
hold on
plot(imag(y),'g*');
hold off;
grid;
title('Tap Value: red: Real channel, green: Imaginary Channel');
pause;

for i=5:5
   %plot(iq_data_e(1+(i-1)*400:-1+i*400),'*');
   plot(iq_data_e(2000:4000),'*');
	grid;
	title('Equalized Constellation');
	axis([-4 4 -4 4]);
	pause;
end


%	plot(iq_data_e(iie-200:iie-1),'*');
%	grid;
%	title('Equalized Constellation');
%	axis([-4 4 -4 4]);
%	%pause;


%ll=30;
%for i_plot=1:5;
%add_off=0;
%offset=10;
%	plot(real(iq_data_in(offset+(i_plot-1)*ll-add_off:offset+i_plot*ll-add_off)),'r*');
%	hold on;
%	plot(real( iq_data_e(offset+(i_plot-1)*ll:offset+i_plot*ll)),'g*');
%	hold off;
%	title('After Equalization: i data (red) and demodulated i data (green)');
%	grid
%%pause
%end	% of i_plot


	plot(real(iq_data_in(iie-30:iie-5)),'r*');
	hold on;
	plot(real( iq_data_e(iie-30:iie-5)),'g*');
	hold off;
	title('After Equalization: i data sent (red) and demodulated i data (green)');
	grid
pause

i_end=converge_i+500;
bit_i=real(bit_error);
bit_q=imag(bit_error);

for ii_e=i_end:iie

	if     bit_i(ii_e)==2
			bit_i(ii_e)=1;
	elseif bit_i(ii_e)==3
			bit_i(ii_e)=2;
	elseif bit_i(ii_e)==4;
			bit_i(ii_e)=1;
	elseif bit_i(ii_e)==5
			bit_i(ii_e)=3;
	elseif bit_i(ii_e)==6;
			bit_i(ii_e)=2;
	end

	if     bit_q(ii_e)==2
			bit_q(ii_e)=1;
	elseif bit_q(ii_e)==3
			bit_q(ii_e)=2;
	elseif bit_q(ii_e)==4;
			bit_q(ii_e)=1;
	elseif bit_q(ii_e)==5
			bit_q(ii_e)=3;
	elseif bit_q(ii_e)==6;
			bit_q(ii_e)=2;
	end
end

bit_error_rate=(sum(bit_i(i_end:iie))+sum(bit_i(i_end:iie)))...
               /max([iie-i_end 1])/8; 		% Eight bits per symbol;
'BER Before FEC'
bit_error_rate

M=256.;   % number of amplitude and phase levels
%m=3*log2(M)/(M-1.);
%u=sqrt(m*10^0.513):sqrt(m*10^0.0005):sqrt(m*10^2.5);
%EbNo=10.*log10(u.^2/m);
%Y=erfc(u/sqrt(2));
%EsNo=EbNo+10*log10(38.236*6.952/41.712);  %5.304 MS/s symbol rate for 6-MHz channel
%A=(1.-(1./sqrt(M)));
%Po=A*Y;
%P1=(2*Po-Po.^2)./log2(M);  % Uncoded 64-QAM BER versus the QAM SNR (in dB)
P1=bit_error_rate;
Ps=P1*log2(M);
% Apply R-S T=8 (204,188) coding to 64-QAM BER
C1=102*203;
C2=34*202*203;
C3=17*101*201*203;
C4=17*40*101*201*203;
C5=17*20*29*67*101*199*198;
C6=5*17*29*67*99*101*197*199;
Pr=1.-(1.-Ps).^204-204.*Ps.*(1-Ps).^203;
Prs= Pr-C1*(Ps.^2).*(1.-Ps).^202-C2*(Ps.^3).*(1.-Ps).^201;
Prr=Prs-C3*(Ps.^4).*(1.-Ps).^200-C4*(Ps.^5).*(1.-Ps).^199;
P2=Prr - C5*(Ps.^7).*(1.-Ps).^197;
P3=P2 -C6*(Ps.^8).*(1.-Ps).^196;
%semilogy(EsNo,P1,'-r',EsNo,P3,'-b')
%axis([20 30 1.e-10 0.1])
%xlabel('64-QAM SNR (dB)')
%B=20:1:36;
%set(gca,'xtick',B)   % set the x-axis tick marks from vector B
%ylabel('Uncoded & RS T=8 (204,188) Coded 64-QAM BER')
%table=[EsNo;P3]'
'BER after FEC'
P3