new_psk_performance.m

分析psk调制与解调的函数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                  %bpsk和Qpsk调制信号在高斯信道下的性能分析%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [cer,ber] = new_psk_performance(n,index,SNR)      
%n,信源产生二进制信号的个数；
 % index：modulation index
 %		1---bpsk
 %		2---qpsk
 %SNR，信噪比
                         
%%%%%%%%%%%%%%%%%%%%%
% 信源产生二进制信号%
%%%%%%%%%%%%%%%%%%%%%
source_output = round(rand(1,n));                   %源输出
%display(source_output);
%%%%%%%%%%%%%%%%%%%%%
%信号mpsk调制%
%%%%%%%%%%%%%%%%%%%%%
input_frame = source_output;
f_length = length(input_frame)/index;                  %调制后信号每帧的长度
psk_input_I = zeros(1,f_length);
psk_input_Q = zeros(1,f_length);
% note: Matlab index starts from 1
switch index
case 1,
    BPSK_I = [-1 1];    % refer to Table82 on page21 of IEEE802.11a
    psk_input_I = BPSK_I(input_frame+1);
    output_digit = psk_input_I ;
case 2,
    QPSK_IQ = [-1 1];   % refer to Table83 on page21 of IEEE802.11a
    psk_input_I = QPSK_IQ(input_frame(1:2:end)+1);
    psk_input_Q = QPSK_IQ(input_frame(2:2:end)+1);
    output_digit = psk_input_I + (psk_input_Q)*j;
end;
%display(output_digit);
%接下来是键控调制到模拟信号
f = 1000;                                   %载波频率
t = 0:.01:f_length;                        
q=pi/3;                                     %载波初始相位
carrier = cos(2*pi*f*t+q);                    %载波
t_length = length(t);
input_modu = zeros(1,t_length);             
m = 0;
  for  t=0:.01:f_length
      m = m+1;
       for i=0:1:f_length-1
           if  (i<=t) &(i+1>t)
              input_modu(m) = output_digit(i+1);           %对每个码元采样 
           end;
       end;
  end;
%display(input_modu);
output_modu = input_modu.*carrier;                  %调制后的输出

%subplot(211);
%plot(t,output_modu);
%display(output_modu);
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%通过高斯信道%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%SNR = -15;                            %信噪比  以db为单位
received_sig=awgn(output_modu,SNR,'measured');            %信号加入高斯白噪声
%subplot(212);
%plot(t,received_sig);
%display(received_sig);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%载波频率和相位估计，利用锁相环，最小二乘法%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t1=0:0.01:1;
f1=f;           %已知信号的频率
q1=q;         %已知信号的初始相位

b=length(t1);      %要求t的长度为奇数
x=real(received_sig);

f2=zeros(1,b);         %定义恢复的频率
q2=zeros(1,b);         %定义恢复的相位
f2(1)=f1;
q2(1)=q1;
M=(b-1)/2;
sum1=zeros(b);
sum2=zeros(b);           
for k=1:1:b-1                  %利用最小二乘（锁相环）估计载波频率和相位
    i=0;
    for a=-M:1:M
        i=i+1;
        sum1(k+1,i+1)=sum1(k,i)+a*x(a+M+1)*sin(2*pi*f2(k)*a+q2(k));
        sum2(k+1,i+1)=sum2(k,i)+1/M*x(a+M+1)*sin(2*pi*f2(k)*a+q2(k));
    end;
    f2(k+1)=f2(k)-3/(4*pi*M^3)*sum1(k+1,end);
    q2(k+1)=q2(k)-1/M*sum2(k+1,end);
end;
%display(sum1);
%display(f1);
%display(f2(m));
%display(q1);
%display(q2(m));
t = 0:.01:f_length; 
carrier_guji=cos(2*pi*f2(b)*t+q2(b));

%plot(t,received_sig,t,carrier_guji,'red');
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%接收端检测%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
%采用双路相关检测结构
T = 1;             %bit interval
tt=0:0.01*T:T;
u=length(tt);
cosref(1,:)=carrier_guji;      %载波信号
cosref(2,:)=-carrier_guji;     %载波反相信号

received_sig_I=real(received_sig);
received_sig_Q=imag(received_sig);

for k=1:f_length
    for i=1:2
        for n=1:u
        code_sig_I(k,i,n)=received_sig_I((k-1)*(u-1)+n)*cosref(i,n);   %每个码元周期内接收信号与载波信号相乘
        code_sig_Q(k,i,n)=received_sig_Q((k-1)*(u-1)+n)*cosref(i,n);
        end;
    end;
end;
for  k=1:f_length
    for i=1:2
        code_sum_I(k,i)=sum(code_sig_I(k,i,:));    %对每个码元周期进行积分，也就是求和
        code_sum_Q(k,i)=sum(code_sig_Q(k,i,:));
    end;
    max_code_sum_I(k)=max(code_sum_I(k,:));       %对两路进行判决，因为先验等概，
                                              %发射信号能量相等，所以判决门限为零
                                              %所以判决就等价于求两路中的最大值
    max_code_sum_Q(k)=max(code_sum_Q(k,:));                                          
end;
for k=1:f_length
    for i=1:2
        if code_sum_I(k,i)==max_code_sum_I(k);
           indexsum_I(k)=i;                    %得到每个码元是载波信号还是载波反相信号 
        end;                                      %i=1,载波信号  i=2，载波反相信号
       if code_sum_Q(k,i)==max_code_sum_Q(k);
           indexsum_Q(k)=i;             
       end;
    end;
end;
reference_one=[1 -1];
for k=1:f_length
    demu_sig_I(k)=reference_one(indexsum_I(k));      %得到解调的信号1或者-1
    demu_sig_Q(k)=reference_one(indexsum_Q(k));
    if index==1
        demu_sig_Q(k)=0;
    end;
    demu_sig(k)=demu_sig_I(k)+j*demu_sig_Q(k);
end;
%display(demu_sig);
reference_two=[1,0];
for k=1:f_length
    emit_sig_I(k)=reference_two((-demu_sig_I(k)+3)/2);    %得到发射的原始信号0或者1
    if index==2
    emit_sig_Q(k)=reference_two((-demu_sig_Q(k)+3)/2);
    end;
end;
  
   if index==1                       %如果是BPSK
        for k=1:f_length
           emit_sig(k)=emit_sig_I(k);
        end;
   end;
   if index==2                       %如果是QPSK
       emit_sig=ones(1,f_length*index);
       for y=1:f_length
          emit_sig(2*y-1)=emit_sig_I(y);
          emit_sig(2*y)=emit_sig_Q(y);   
       end;
   end;
%display(emit_sig);
     
for k=1:f_length
    error_code(k)=demu_sig(k)-output_digit(k);       %比较解调的码元与调制得到的码元
end;
for k=1:f_length*index
    error_bit(k)=emit_sig(k)-source_output(k);       %比较解调的比特信号与原始信号
end;

num1=0;
num2=0;
for k=1:f_length
    if error_code(k)==0
        num1=num1+1;                    %正确解调的码元数
    end;
end;
for k=1:f_length*index
    if error_bit(k)==0
        num2=num2+1;                    %正确解调的比特数
    end;
end;
num_error_code=f_length-num1;          %错误码元数
num_error_bit=f_length*index-num2;           %错误比特数
format short e;
cer=num_error_code/f_length;           %误码率
ber=num_error_bit/(f_length*index);            %误比特率
%display(cer);
%display(ber);
format;