msk_noise_c_dj.m

本程序为通信系统中的MSK调制程序

%MSK调制，解调
clear all
close all

%--------------------------------------------------------------------------
%参数设置
data_len = 30000;           %码元个数
sample_number = 8;          %采样个数
Rb = 24000;                 %码元速率
fc = 96000;                 %载波频率
%**************************************************************************

%--------------------------------------------------------------------------
%随机产生传输信号
data = rand_binary(data_len);
%**************************************************************************

%--------------------------------------------------------------------------
%MSK基带调制
[signal_out,I_out,Q_out] = mod_msk(data,data_len,sample_number,Rb);
%**************************************************************************

%--------------------------------------------------------------------------
%中频搬移
multi = fc/Rb;
I_temp=interp(I_out,multi);
Q_temp=interp(Q_out,multi);

Fs=fc*sample_number;
t=1/Fs:1/Fs:length(I_temp)*1/Fs;
signal_i=I_temp.*cos(2*pi*fc*t);
signal_q=Q_temp.*sin(2*pi*fc*t);
signal_mod=I_temp.*cos(2*pi*fc*t)-Q_temp.*sin(2*pi*fc*t);
% signal_mod1=real((I_temp+j*Q_temp).*exp(j*2*pi*fc*t));
%**************************************************************************

amp=[1 0 0.69 0 0.24];
dj=length(amp);
for i = 1:dj
    signal_mod=signal_mod+amp(i)*delay(signal_mod,i-1,sample_number);
end

%--------------------------------------------------------------------------
%加噪声
for SNR = 0:10
% signal_mod1 = awgn(signal_mod,SNR);

noise=randn(1,length(signal_mod));
b=5/50;c=6/50;d=19/50;e=20/50;
F1=[0,b,c,d,e,1];
A1=[0,0,1,1,0,0];
bpf=firls(300,F1,A1);

noise1=conv(noise,bpf);
noise_out=noise1(150+1:150+length(signal_mod));

Power_noise=sum(noise_out.^2);
Power_signal=sum(signal_mod.^2);
Coe_SNR=sqrt(1*Power_signal/Power_noise/(10^(SNR/10)));

signal_mod1 = signal_mod+(Coe_SNR)*noise_out;    %加入噪声以后的信号

    %--------------------------------------------------------------------------
    %去载波
    N=300;                                              % 滤波器的阶数为(N+1)  
    F=[0,fc-1000,fc+1000,Fs/2]*2/Fs;
    A=[1,1,0,0];
    lpf=firls(N,F,A);
    [amp_lpf,w]=freqz(lpf);
    
    I_dem=signal_mod1.*cos(2*pi*fc*t)*2 .* cos(pi*t*Rb/2);
    I_dem=conv(I_dem,lpf);
    I_dem=I_dem(N/2+1:N/2+length(I_temp));
    Q_dem=signal_mod1.*sin(2*pi*fc*t)*2 .* sin(pi*t*Rb/2);
    Q_dem=conv(Q_dem,lpf);
    Q_dem=-Q_dem(N/2+1:N/2+length(I_temp));
    
    I_dem_out=zeros(1,length(I_dem)/multi);         % 抽取
    Q_dem_out=zeros(1,length(Q_dem)/multi);
    for i=1:length(I_dem_out)
       I_dem_out(i)=I_dem(multi*(i-1)+1);
       Q_dem_out(i)=Q_dem(multi*(i-1)+1);
    end;    
    %**************************************************************************

    %--------------------------------------------------------------------------
    %解调，判决
    demod_data = zeros(1,data_len);
    for i = 1:data_len
        I_dem1(i) = I_dem_out(i*sample_number)>0;
        Q_dem1(i) = Q_dem_out(i*sample_number)>0;
    end
    for i = 1:2:data_len
        m1(i) = Q_dem1(i);
    end
    for i = 2:2:data_len
        m1(i) = I_dem1(i);
    end
    m1 = 2*m1-1;
    %**************************************************************************

    %--------------------------------------------------------------------------
    %差分解码
    demod_data = zeros(1,data_len);
    demod_data(1) = 1 * m1(1);
    for i = 2:data_len
        demod_data(i) = m1(i-1) * m1(i);
    end
    %**************************************************************************
    
    %--------------------------------------------------------------------------
    %计算误码率
    [num,ber(SNR+1)]=symerr(demod_data,data);
    %**************************************************************************
end              % /for snr
%**************************************************************************

%--------------------------------------------------------------------------
%误码率曲线
semilogy([0:10],ber,'r*-');
%**************************************************************************

%--------------------------------------------------------------------------
%误码率理论值
snr = 0:0.1:10;
for i = 1:length(snr)
    snr1(1,i) = 10^(snr(1,i)/10);
    ps(1,i) = 1/2 * erfc(sqrt(snr1(1,i)));
    pe(1,i) = 2 * ps(1,i);
end
hold on
semilogy([0:.1:10],pe);
%**************************************************************************