qpsk_rayleigh.m

qpsk在Rayleigh下的调制和解调

% qpsk_rayleigh.m
%
% Simulation program to realize QPSK transmission system over Rayleigh channel
%

%******************** Preparation part **************************************

clear;
clc;

sr=256000.0; % Symbol rate
ml=2;        % ml:Number of modulation levels (BPSK:ml=1, QPSK:ml=2, 16QAM:ml=4)
br=sr.*ml; % Bit rate
nd=100;   % Number of symbols that simulates in each loop
ebn0=10;      % Eb/N0
IPOINT=8;    % Number of oversamples
sample2=32;

%************************* Filter initialization ***************************

irfn=21;                  % Number of taps
alfs=0.5;                 % Rolloff factor
[xh] = hrollfcoef(irfn,IPOINT,sr,alfs,1);   %Transmitter filter coefficients 
[xh2] = hrollfcoef(irfn,IPOINT,sr,alfs,0);  %Receiver filter coefficients 

%******************** START CALCULATION *************************************

nloop=100;  % Max number of simulation loops
noe = 0;    % Number of error data
nod = 0;    % Number of transmitted data
j=sqrt(-1);
for iii=1:nloop
    
%*************************** Data generation ********************************  
	
    data1=rand(1,nd*ml)>0.5;  % rand: built in function
    
%*************************** QPSK Modulation ********************************  

    [ich,qch]=qpskmod(data1,1,nd,ml);
	[ich1,qch1]= compoversamp(ich,qch,length(ich),IPOINT); 
	[ich2,qch2]= compconv(ich1,qch1,xh); 
    R1=ich2+j.*qch2;
    Mag1=abs(R1);
%****************************Add signals to carriers*********************

ich_sample2=oversample2(ich2,sample2);
qch_sample2=oversample2(qch2,sample2);

%Frequency of carriers is 4 times as large as 
%that of singnals after oversample(IPOINT=8)
cos_carrier=cos((8*pi/sample2)*[1:length(ich_sample2)]);
sin_carrier=sin((8*pi/sample2)*[1:length(qch_sample2)]);
                   
ich_sample2=ich_sample2.*cos_carrier;
qch_sample2=qch_sample2.*sin_carrier;

%*************************** Fading Channel ******************************

    %  define a frequency vector and a magnitude vector to simulate the classic 'bathtub' shape
  
    f=(0.0:0.05:1.0);
    m=[1.0,1.2,1.5,1.9,2.8,4.0,6.0,9.0,15.0,25.0,1.00,0.05,0.03,0.02,0.01,0.005,0.005,0.005,0.005,0.005,0.005];
  
    B=fir2(16,f,m);		% design an FIR filter based on the shape above. 
    ich_sample2_R=filter(B,1,ich_sample2);
    qch_sample2_R=filter(B,1,qch_sample2);
   
%**************************** Attenuation Calculation ***********************
	
    spow=sum(ich2.*ich2+qch2.*qch2)/nd;  % sum: built in function
	attn=0.5*spow*sr/br*10.^(-ebn0/10);
	attn=sqrt(attn);  % sqrt: built in function
     
%********************* Add White Gaussian Noise (AWGN) **********************
	
    [ich31,qch31]= comb(ich_sample2_R,qch_sample2_R,attn);% add white gaussian noise
     R2=ich31+j.*qch31;
     Mag2=abs(R2);
%********************Resume signals from carriers**************************

    ich3=ich31.*cos_carrier;
    qch3=qch31.*sin_carrier;
    
    %butterworth filter(cut-off frequency is chosen as same as carrier)
    [b,a]=butter(10,0.25);
    
    %filter the signals
    ich3=filter(b,a,ich3);
    qch3=filter(b,a,qch3);
    
    %resume dada from oversampled signals every 32 points
    ich3_data=ich3(1:sample2:length(ich3));
    qch3_data=qch3(1:sample2:length(qch3));
    
    [ich4,qch4]= compconv(ich3_data,qch3_data,xh2);

    syncpoint=irfn*IPOINT+1;
    ich5=ich4(syncpoint:IPOINT:length(ich4));
    qch5=qch4(syncpoint:IPOINT:length(qch4));
        
%**************************** QPSK Demodulation *****************************

    [demodata]=qpskdemod(ich5,qch5,1,nd,ml);

%************************** Bit Error Rate (BER) ****************************

  noe2=sum(abs(data1-demodata));  % sum: built in function
  nod2=length(data1);  % length: built in function
  noe=noe+noe2;
  nod=nod+nod2;
fprintf('%d\t%e\n',iii,noe2/nod2);
%************************Check if noe>=100******************************

 % if noe>=10
  %   flag=1;
   %  break;
  %end
%fprintf('%d\t%e\n',iii,noe2/nod2);  % fprintf: built in function

end  %for iii=1:nloop

%********************** Output result ***************************
%if flag==1
   ber = noe/nod;
   fprintf('%d\t%d\t%d\t%e\n',ebn0,noe,nod,noe/nod);  % fprintf: built in function
%else
  % fprintf('Time out!\n');
%end
 %end
%*********************Plot the data of the last loop**************

figure(1);

subplot(2,2,1);
plot(data1(1:10),'c');
title('调制信号');
subplot(2,2,2);
DATA1=abs(fft(data1));
plot(fftshift(DATA1),'r');
title('调制信号频谱');

subplot(2,2,3);
plot(ich_sample2,'b');
title('已调信号(信道I)');
subplot(2,2,4);
ICH_SAMPLE2=abs(fft(ich_sample2));
plot(fftshift(ICH_SAMPLE2),'y');
title('已调信号频谱(信道I) ');

figure(2);

subplot(2,2,1);
plot(ich31,'k');
title('接收信号(信道I)');
subplot(2,2,2);
ICH31=abs(fft(ich31));
plot(fftshift(ICH31),'g');
title('接收信号频谱(信道I)');

subplot(2,2,3);
plot(demodata(1:10),'g');
title('解调信号(前10周期)');
subplot(2,2,4);
DEMODATA=abs(fft(demodata));
plot(fftshift(DEMODATA),'r');
title('解调信号频谱');
figure(3);
subplot(2,1,1);
plot(Mag1);
subplot(2,1,2);
plot(fftshift(Mag1));

figure(4);
plot(Mag2);

scatterplot([ich' qch'],'b');
title('传输前信号的星座图');

scatterplot([ich5' qch5'],'m');
title('S传输信号的星座图');

%******************** end of file ***************************