qam256.m

本文件内包含了QPSKandQAM256的系统仿真

%*************************************************
%      256 QAM EQUALIZER TESTING FILE             *
%************************************************* 
%						 *
%      	DATE:   5/28/98; 			 *
%       AUTHOR: QIN ZHANG 	 		 *
%         					 *	 
%	FILE NAME: EQUA_256.M  (12-BIT 256 QAM)	 *
%	C:QAM\SIMULATION                         *  
%  transmitter alpha=.12 receiver alpha=.12  *
%************************************************* 
%
%	Two Channel Third Order Simulation;
%path(path,'c:\matlab\mat_apl');   
%path(path,'c:\matlab\mg');  
%path(path,'c:\matlab\data');  
clear;

fs =5e6*4; %5.304e6*4;   			 %42.884296e6/2;		  % 20.227640MHz sampling frequency;
f_symbol=fs/4;				  % 10.113820MSPS

taps_t=5*2*2*2*2*2;		% 160-tap filter;
alpha_t=0.1865;			%alpha_t=0.1865	% Transmiter SRRC filter alpha;

taps_r=taps_t; 			% receiver filter taps;
alpha_r=alpha_t;			% ideal: 0.1865	% receiver fitler alpha;

n_filter=taps_t*2-1; 	% raised cosine filter taps;
n_bit_rc =16;				% raised cosine filter coefficient bit;
n=1024*4*2;					% random data size;
n_demo_filter=6;			% demodulation low pass filter taps;
f_cutoff_demo=5e6;		% demodulation low pass cutoff frequency;
n_bit_da=16;				% D/A converter resolution;
n_bit_sin=16;				% Sine and Cosine Coefficient;
ph_error=0;             % carrier phase error;
ny_full=1;					% full nyquist filter; 0 partial
ny_n=10;           		% partial nyquist filter, 8 taps on left and right from center tap;

%*****************************************************************
%    1)  Generate Square Root Raised Cosine Filter Coefficients  *
%*****************************************************************


t=1:taps_t;
t=(t-floor(taps_t/2)-1+1e-10)/fs;

hrc_t=(sin(pi*t*f_symbol*(1-alpha_t)) + 4*alpha_t*t*f_symbol.*...
       cos(pi*t*f_symbol*(1+alpha_t))    )./...
       ( pi*t*f_symbol.*(1-(4*alpha_t*t*f_symbol).^2));

%hrc_t=hrc_t.*blackman(taps_t)';
hrc_t=hrc_t.*kaiser(taps_t,6)';

plot(t, hrc_t);
grid;
title('hrc_t impulse response');
pause;

hrc_r=hrc_t;
hrc_t=hrc_t/sum(hrc_t);
hrc_r=hrc_r/sum(hrc_r);

maxmin=max(abs(hrc_t));
module=1/(2^(n_bit_rc-1)-1)*maxmin; 	% define LSB;
hrc_t_i=round(hrc_t/module);

% transform the fitler to polyphase;

hrc_t_4=reshape(hrc_t,4,taps_t/4);
hrc_t_4(3,:)=-hrc_t_4(3,:);
hrc_t_4(4,:)=-hrc_t_4(4,:);

%****************************************
%      3) Get the random data		*	
%****************************************
total_fft=20;

for fft_loop=1:total_fft;
   
for rand_data_set=1:2

rand('state', sum(100*clock));
[idat,qdat] = Pn_iq4(n,1);	% -7.5 to 7.5;
%idat=ones(size(idat));
%qdat=idat;
	
	plot(idat(1:20),'*');
	grid;
	title('I Data');
%pause;
	
	iq_data=idat+qdat*1i;
	plot(iq_data(100:900),'*');
	grid;
	title('I/Q Data Constellation');
	axis([-4 4 -4 4]*2);
	pause
	% create a test impulse
	%idat=zeros(size(idat));
	%idat(200)=8;
	%qdat(203)=-8;
	%idat(199)=-8;

%*************************************************
%        4) Base Band Filtering & modulation     *
%*************************************************

%grp_delay=taps_t/2;
offset=801+800;		
nnnn=1024*4;
	data_fft=10*log10(abs(fftshift(fft(kaiser(1024*4,11.)'.*...
       (idat(offset:offset+1024*4-1)+qdat(offset:offset+1024*4-...
 	1))))));
%	data_fft=20*log10(abs(fftshift(fft(kaiser(nnnn,11.)'.*...
%	idat(offset:offset+nnnn-1)))));
	data_fft=data_fft-max(data_fft);
	ff=-nnnn/2:nnnn/2-1;
	ff=ff/nnnn*2*fs/2;
       plot(ff,data_fft,'b');
	title(' I/Q Data Spectrum Before Filtering ');
	axis([0, 1e7 -85 0]);
	grid;
	xlabel('FREQUENCY (HZ)');
	ylabel('MAGNITUDE (dB) ');
pause;
%	print fig2 -dmeta;

%
%    polyphase filtering with Post Filtering Carrier Recombination;
%

if_data1=4*filter(hrc_t_4(1,:),1,idat+0.5)-4*sum(hrc_t_4(1,:))*.5;
if_data2=4*filter(hrc_t_4(2,:),1,qdat+0.5)-4*sum(hrc_t_4(2,:))*.5;
if_data3=4*filter(hrc_t_4(3,:),1,idat+0.5)-4*sum(hrc_t_4(3,:))*.5;
if_data4=4*filter(hrc_t_4(4,:),1,qdat+0.5)-4*sum(hrc_t_4(4,:))*.5;

if_data=[if_data1; if_data2; if_data3; if_data4];
rf_data=reshape(if_data,1,n*4);

offset=200;		 	% filter warm up time;

%	 plot spectrum;
	nnnn= 1024*8;
	data_fft=20*log10(abs(fftshift(fft(kaiser(nnnn,11.)'.*...
				 rf_data(offset:offset+nnnn-1)))));
	data_fft=data_fft-max(data_fft);
	ff=-nnnn/2:nnnn/2-1;
	ff=ff/nnnn*2*fs/2;
  plot(ff,data_fft);
   title(' I/Q Data Spectrum After Filtering ');
	axis([0, 1e7 -100 0]);
	grid;
	xlabel('FREQUENCY (HZ)');
	ylabel('MAGNITUDE (dB) ');
	pause;
%	print fig3 -dmeta;
%       axis([3e6 3.5e6 -90 0]);
%	%%pause;



%************************************************
%	Interpolation				fs=40	MHz				*
%************************************************

rf_data0=zeros(size(rf_data));
rf_data_8=[rf_data; rf_data0];
rf_data_8=reshape(rf_data_8,1,n*4*2);

fs_8=fs*2;

n_filter=160;
f=[ 0 .41   .56             1];
m=[1   1   10^(-6) 10^(-8) ];

b=fir2(n_filter, f, m);
%[H,F] = FREQZ(b,1, 1024,fs_8);
%plot(F, 20*log10(abs(H)));
%grid;
%pause

rf_data_8=filter(b,1,rf_data_8);

%	 plot spectrum;
%	nnnn= 1024*8;
%	data_fft=20*log10(abs(fftshift(fft(kaiser(nnnn,11.)'.*...
%				 rf_data_8(offset:offset+nnnn-1)))));
%  data_fft=data_fft-max(data_fft);
%	ff=-nnnn/2:nnnn/2-1;
%	ff=ff/nnnn*2*fs_8/2;
%	plot(ff,data_fft);
%	title(' I/Q Data Spectrum After Interplation x2');
%	axis([0 max(ff) -100 0]);
%	grid;
%	xlabel('FREQUENCY (HZ)');
%	ylabel('MAGNITUDE (dB) ');
%	pause;



%************************************************
%	Interpolation	fs=80 MHz							*
%************************************************

rf_data0=zeros(size(rf_data_8));
rf_data_16=[rf_data_8; rf_data0];
rf_data_16=reshape(rf_data_16,1,n*4*2*2);

fs_16=fs*2*2;

n_filter=100;
f=[ 0 .22   .77             1];
m=[1   1   10^(-10) 10^(-10) ];

b=fir2(n_filter, f, m);
%[H,F] = FREQZ(b,1, 1024,fs_16);
%plot(F, 20*log10(abs(H)));
%grid;
%pause

rf_data_16=filter(b,1,rf_data_16);

%	 plot spectrum;
%	nnnn= 1024*8;
%	data_fft=20*log10(abs(fftshift(fft(kaiser(nnnn,11.)'.*...
%				 rf_data_16(offset:offset+nnnn-1)))));
%   data_fft=data_fft-max(data_fft);
%	ff=-nnnn/2:nnnn/2-1;
%	ff=ff/nnnn*2*fs_16/2;
%   plot(ff,data_fft);
%   title(' I/Q Data Spectrum After Interplation X4');
%	axis([0 max(ff) -100 0]);
%	grid;
%	xlabel('FREQUENCY (HZ)');
%	ylabel('MAGNITUDE (dB) ');
 %  pause;
%************************************************
%	Interpolation	fs=160 MHz							*
%************************************************

rf_data0=zeros(size(rf_data_16));
rf_data_32=[rf_data_16; rf_data0];
rf_data_32=reshape(rf_data_32,1,n*4*2*2*2);

fs_32=fs*2*2*2;

n_filter=60;
f=[ 0 .11   .77             1];
m=[1   1   10^(-10) 10^(-10) ];

b=fir2(n_filter, f, m);
%[H,F] = FREQZ(b,1, 1024,fs_32);
%plot(F, 20*log10(abs(H)));
%grid;
%pause

rf_data_32=filter(b,1,rf_data_32);

%	 plot spectrum;
%	nnnn= 1024*8;
%	data_fft=20*log10(abs(fftshift(fft(kaiser(nnnn,11.)'.*...
%				 rf_data_32(offset:offset+nnnn-1)))));
%	data_fft=data_fft-max(data_fft);
%	ff=-nnnn/2:nnnn/2-1;
%	ff=ff/nnnn*2*fs_32/2;
%	plot(ff,data_fft);
%	title(' I/Q Data Spectrum After Interplation X8');
%	axis([0 max(ff) -100 0]);
%	grid;
%	xlabel('FREQUENCY (HZ)');
%	ylabel('MAGNITUDE (dB) ');
%	pause;
   
if rand_data_set == 1   
%************************************************
%	UP converting to 17.25MHz							*
%************************************************   
   
t=[0:n*32-1]/fs_32;
c_cos=cos(2*pi*(27.25-5-.2)*1e6*t); 	% carrier sine;

rf_data1=rf_data_32.*c_cos;	% up converting


n_filter=351;
f=[ 0 2/8   2.2/8     3/8  1];
m=[1   1   10^(-10) 10^(-10) .001 ];

b=fir2(n_filter, f, m);
%[H,F] = FREQZ(b,1, 1024,fs_32);
%plot(F, 20*log10(abs(H)));
%grid;
%pause

rf_data1=filter(b,1,rf_data1);

%	 plot spectrum;
%	nnnn= 1024*8;
%	data_fft=20*log10(abs(fftshift(fft(kaiser(nnnn,11.)'.*...
%				 rf_data1(offset:offset+nnnn-1)))));
%   data_fft=data_fft-max(data_fft);
%	ff=-nnnn/2:nnnn/2-1;
%	ff=ff/nnnn*2*fs_32/2;
%   plot(ff,data_fft);
%   title(' RF Data Spectrum After Up Converting');
%	axis([0 max(ff) -100 0]);
%	grid;
%	xlabel('FREQUENCY (HZ)');
%	ylabel('MAGNITUDE (dB) ');
%  pause;
   
else
   
%************************************************
%	UP converting to 20.25MHz							*
%************************************************   
   
t=[0:n*32-1]/fs_32;
c_cos=cos(2*pi*(27.25-5+5)*1e6*t); 	% carrier sine;

rf_data2=rf_data_32.*c_cos;	% up converting


n_filter=350;
f=[ 0 2.5/8   2.7/8             1];
m=[1   1   10^(-10) 10^(-10) ];

b=fir2(n_filter, f, m);
%[H,F] = FREQZ(b,1, 1024,fs_32);
%plot(F, 20*log10(abs(H)));
%grid;
%pause

rf_data2=filter(b,1,rf_data2);

%	 plot spectrum;
%	nnnn= 1024*8;
%	data_fft=20*log10(abs(fftshift(fft(kaiser(nnnn,11.)'.*...
%				 rf_data2(offset:offset+nnnn-1)))));
%   data_fft=data_fft-max(data_fft);
%	ff=-nnnn/2:nnnn/2-1;
%	ff=ff/nnnn*2*fs_32/2;
%   plot(ff,data_fft);
%   title(' RF Data2 Spectrum After Up Converting');
%	axis([0 max(ff) -100 0]);
%	grid;
%	xlabel('FREQUENCY (HZ)');
%	ylabel('MAGNITUDE (dB) ');
%	pause;
end

end % of rand_data_set

rf_two_channel(fft_loop,:)=(rf_data1+rf_data2);
%rf_two_channel=rf_data2;

%	 plot spectrum;
	nnnn= 1024*2*2*2*2;
%	data_fft(fft_loop,:)=abs(fftshift(fft(kaiser(nnnn,11.)'.*...
%      rf_two_channel(fft_loop, offset:offset+nnnn-1))));

end % end of fft_loop;

%temp=data_fft(1,:);
%for i=2:fft_loop
%   temp=temp+data_fft(i,:);
%end
%temp=temp/fft_loop;

%   data_fft_db=20*log10(temp);
%   data_fft=data_fft-max(data_fft);
%	ff=-nnnn/2:nnnn/2-1;
%	ff=ff/nnnn*2*fs_32/2;
%   plot(ff/1e6,data_fft_db);
%   title(' Two-Channel 256QAM RF  Spectrum');
%	axis([0 max(ff)/1e6/2 -40 50]);
%	grid;
%	xlabel('FREQUENCY (MHZ)');
%   ylabel('AMPLITUDE (dB) ');
   %print fig_two_tone_QAM -dmeta;

%   pause;
   amp=2;
   for fft_loop=1:total_fft
   ctb(fft_loop,:)=0.01*((amp*rf_two_channel(fft_loop,:)).^3)+amp*rf_two_channel(fft_loop,:);
   
   data_fft(fft_loop,:)=abs(fftshift(fft(kaiser(nnnn,11.)'.*...
      ctb(fft_loop, offset:offset+nnnn-1))));
end % end of fft_loop;

temp=data_fft(1,:);
for i=2:total_fft
   temp=temp+data_fft(i,:);
end
temp=temp/total_fft;

   data_fft_db=20*log10(temp);

   ff=-nnnn/2:nnnn/2-1;
	ff=ff/nnnn*2*fs_32/2;
   plot(ff/1e6,data_fft_db, 'r');
   title(' Third Order Distortion');
	axis([0 max(ff)/1e6/2 -30 50]);
	grid;
	xlabel('FREQUENCY (MHZ)');
	ylabel('MAGNITUDE (dB) ');
   pause;
      %print fig_3rd_QAM_2 -dmeta;









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

%************************************************
% Define Channel Model				*
%************************************************


ch_delay=6-1;  	% 6= channel length *4 from the center tap

channel=[ -0.0 0. 0. 0. -0. 0. 0. 0. -0.0002*8*1.1 0. 0. 0. -0.0 0. 0. 0. ...  % pre echo,tap index: [19  18  17  16] 
         0.0 0. 0. 0. -0.0 0. -0.0 0. ...					% tap index: [15 14] 
	1 ...							% center tap	     [13]
	  0.0 0.  0.0 -0.1 0. 0. 0. -0.0 0.0 0. 0. -0. 0. 0. 0. -0. ... %   index: [12 11 10  9] 
          0. 0. 0.  -0.0 0. 0. 0. -0. 0.0 0. 0. -0. 0. 0. 0. -0. ...   % tap index: [8  7  6  5] 
          0. 0. 0.  -0.+i*.0005120 0. 0. 0. -0. 0.0 0. 0. -0. 0. 0. 0. -0. ...   % tap index: [4  3  2  1] 
	  0. 0. 0.0  -0.0 0. 0. 0. -0. 0.0 0. 0. -0. 0. 0. 0. -0. ];
training=0;	% if training=1 use training.
converge_i=0;
test_tap1=13-4;       %  start tap index of the equalizer.  The index is defined in the above  table [1 to 19 ]
test_tap2=13+4;	   %  end   tap index of the equalizer.
i_module=100;
end_iie=1024*4.5;	% end of iie;
mu=.0015*5/2; %1/2^6;	%0.03
n_test=12; 	% Feed Backward taps;
rf_data=filter(channel,1,rf_data);
%****************************************
%		Set the SNR  		*
%****************************************
n_s=.2*2*1.2*(rand(size(rf_data))-.5); 

%*************************************************
%        5)  Demodulation                        *
%*************************************************


t=[0:n*4-1]/fs;
c_sin=sin(2*pi*fs*(1+.0e-3/5)/4*t+0/180*pi); 	% carrier sine;
c_cos=cos(2*pi*fs*(1+.0e-3/5)/4*t+0/180*pi); 	% carrier cosine;

rf_i_d=rf_data.*c_cos;	% down converting
rf_q_d=rf_data.*c_sin;	% down converting

noise_i=n_s.*c_cos;	% down converting
noise_q=n_s.*c_sin;	% down converting


%	% plot spectrum;
%	data_fft=20*log10(abs(fftshift(fft(kaiser(1024,11.)'.*...
%		 rf_i_d(offset:offset+1024-1)))));