awgn_cdma_pe加频差估计.m

这是个扩频系统中的PN码仿真程序

%通信信号处理作业
%
%
%
%
clear;clc;%close all;采用了估计频偏的方法
SNR=10;
% function pe_awgn=awng_cdma_pe(SNR)
err_num=0;
data_bit_length=0;
while err_num<20
    N=200;
    Lc=15;
    Np=8;
    fd=5.6;
    rate=2400;
    if_frequency=1800;%数字上变频
    fs=14400;
    j=sqrt(-1);
    d=1;
    rand('state',sum(100*clock));
    data_bit=round(rand(1,N));
    modula_data=QPSK_Modulate(data_bit);
    modua_data_I=real(modula_data);
    modua_data_Q=imag(modula_data);
    
    code_type=2;  %0 = random 1 = hadamard 2 = gold生成扩品吗
    code_seed=rand('state'); 
    codes=(cdma_codes(1,0,code_type,Lc,code_seed))'; %cdma_codes(K_in,K_out,code_type,L,code_seed) 
    % codes=[1 1 0 0 0 0 1 0 1 1 0 1 1 1 0];
    I_spread=kron(modua_data_I,codes);%产生两路扩频信号
    Q_spread=kron(modua_data_Q,codes);
    over_sample=[1 0 0 0 0 0];%过抽样，消除isi
    I_spread_samp=kron(I_spread,over_sample);
    Q_spread_samp=kron(Q_spread,over_sample);
    num=rcosfir(0.2,4,6,1/rate,'sqrt');%成型滤波
    I_rcosed_delay=conv(I_spread_samp,num);
    Q_rcosed_delay=conv(Q_spread_samp,num);
    
%     for i=1:length(I_spread_samp)
%         I_rcosed(i)=I_rcosed_delay(i+36);
%         Q_rcosed(i)=Q_rcosed_delay(i+36);
%     end

    i=1:length(I_rcosed_delay);
    IF_fre_I=cos(2*pi*if_frequency*(i-1)*(1/fs));%上变频 data_f=14.4/522,所以1.8*522/14.4＝65.25就是中心载波所在的位置
    IF_fre_Q=sin(2*pi*if_frequency*(i-1)*(1/fs));
    I_if=I_rcosed_delay.*IF_fre_I;
    Q_if=Q_rcosed_delay.*IF_fre_Q;
    S_I_Q=I_if+Q_if;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %信道
    %加入噪声
%     SNR=10;
%     E=1;%归一化符号号能量
%     E_chip=1;
%     E_bit=E/2;
%     snr=10^(SNR/10);
%     sigma=sqrt(E_chip/(2*snr*4/15));%noise variance  (snr_in bit=E_bit/(N0))
%     randn('state',sum(100*clock)+50);
%     Noise=sigma*randn(size(S_I_Q));
%     S_I_Q=S_I_Q+Noise;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % %数字下变频
%     IF_fre_I=cos(2*pi*(if_frequency+1.6)*(i-1)*(1/fs));
%     IF_fre_Q=sin(2*pi*(if_frequency+1.6)*(i-1)*(1/fs));

    I_down=S_I_Q.*IF_fre_I;
%     f1=2*rate/fs;
%     low_filter=fir1(34,f1);
%     after_low_tem=conv(I_down,low_filter);
%     for i=1:length(I_down)
%         after_low_I(i)=after_low_tem(i+17);
%     end
    Q_down=S_I_Q.*IF_fre_Q;%直接做就可以了
%     after_low_tem=conv(Q_down,low_filter);
%     for i=1:length(Q_down)
%         after_low_Q(i)=after_low_tem(i+17);
%     end
    
    I_gr=conv((I_down),num);
    Q_gr=conv((Q_down),num);
    for i=1:length(I_spread_samp)
        I_after_gr(i)=I_gr(i+48);%周期为48和36这是为什么？
        Q_after_gr(i)=Q_gr(i+36);
    end
    I_down_samp=downsample(I_after_gr,6);%DOWNSAMPLE 函数没有给出
    Q_down_samp=downsample(Q_after_gr,6);
    
%     i=1;%和上两句完成的功能一样
%     while (i*6)<=length(I_spread_samp)
%         I_down_samp(i)=I_after_gr((i-1)*6+1);
%         Q_down_samp(i)=Q_after_gr((i-1)*6+1);
%         i=i+1;
%     end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算该点的SNR
    
% %     SNR=10;如果使用统计信号的能量，可以分别计算每支路的方差
% 
%     E=(sum(I_down_samp.^2))/length(I_down_samp);%认为能量没有衰减
%     snr=10^(SNR/10);
%     sigma=sqrt(E/(2*snr));%用的是比特信噪比
%     sigma_I=sigma;%正交支路方差为总的一半
%     sigma_Q=sigma;
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
   %下抽之后chip信噪比，转换成比特是为了比较方便
    E=1;%(sum(I_down_samp.^2))/length(I_down_samp);认为能量没有衰减
    snr=10^(SNR/10);
    sigma=sqrt(E/(4*snr));%用的是比特信噪比
    sigma_I=sigma/2;%正交支路方差为总的一半
    sigma_Q=sigma/2;
%     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    randn('state',sum(100*clock)+50);
    Noise_I=sigma_I*randn(size(I_down_samp));
    Noise_Q=sigma_Q*randn(size(Q_down_samp));
    I_down_samp=I_down_samp+Noise_I;
    Q_down_samp=Q_down_samp+Noise_Q;

    % %验证
    % aa=fft(I_rcosed);
    % dd=fft(Q_rcosed);
    % bb=fft(I_down);
    % cc=fft(Q_down);
    % figure(1);
    % plot(abs(aa));
    % figure(2);
    % plot(abs(bb));
    % figure(3);
    % plot(abs(cc));
    % figure(4);
    % plot((abs(dd)));
    data_length=length(I_down_samp);
    length_symbol=round(data_length/Lc);
    for i=1:length_symbol
        tem=I_down_samp((i-1)*(Lc)+1:(i-1)*(Lc)+15);
        tem_sum=sum(tem.*codes);
        symbol_I(i)=tem_sum*2/Lc;
        tem=Q_down_samp((i-1)*(Lc)+1:(i-1)*(Lc)+15);
        tem_sum=sum(tem.*codes);
        symbol_Q(i)=tem_sum*2/Lc;
    end
    symbol=symbol_I+j*symbol_Q;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %相关解调，考虑了频偏估计，结果证明，加估计比不加性能要hao，
    Zn=symbol;
    pliot=modula_data(1:Np);
    symbol=Mpsk_symbol_dec(pliot,Zn,Np);
    symbol_demoulate=qpsk_bit_mapping(symbol,d);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     symbol_demoulate=QPSK_Demodulate(symbol,d);%%%%不加估计时候使用
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    err_num1=pe_bit(data_bit,symbol_demoulate);
    err_num=err_num+err_num1
    data_bit_length=data_bit_length+length(data_bit);
end
pe_awgn=err_num/data_bit_length;
% pe_awgn=pe_bit(data_bit,symbol_demoulate);