lfh_dpsk.m

慢速跳频

clc,clear
close all

% 慢跳频系统
% 数据采用 DPSK 调制


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%         跳频参数设置             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ts = 0.001;          % 数据比特周期 (s)

fm = 50000;          % DPSK 信号载频中心频率
fs = 1/Ts;           % 数据传输速率
Bs = 2*fs;           % DPSK 信号带宽 (跳频频率最小间隔)

Tc = 4*Ts;           % 跳频周期 (Tc>=Ts)
K = 16;              % 跳频频率个数
fh = 100000;         % 跳频信号中心频率
Bt = K*Bs;          % 发射信号带宽

t = 0.1;              % 仿真时间
fsample = 800000;   % 采样速率 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%         跳频信号产生             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Data_len = floor(t/Ts);                                     % 用户数据长度
Data_diff = randsrc(Data_len,1,(0:1));                      % 用户数据的差分数据
Data = xor(Data_diff(1:end-1),Data_diff(2:end));            % 用户数据 
%%%%%%%%%%%%%%%%        MFSK 调制       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
St_baseband = DMOD(Data_diff,fm, fs, fsample,'psk', 2);         % PSK 调制信号
% figure, plot(St_baseband),title('DPSK基带信号');
signal_len = length(St_baseband);                           % 仿真时间内采样数据长度
%%%%%%%%%%%%%%%%%%%   跳频载波产生   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FH_seq = randsrc(ceil(signal_len/fsample/Tc),1,(0:K-1));         % 跳频序列
FH_sin = zeros(signal_len,1);                               % 本地跳频信号载波
phase = 2*pi*rand;
for i = 1:ceil(signal_len/fsample/Tc)
    f2 = fh + (FH_seq(i)*2 - (K-1))*Bs/2;                   % 跳频频点中心频率
    FH_sin((i-1)*fsample*Tc + 1:i*fsample*Tc) = cos(2*pi*f2*(0:fsample*Tc-1)/fsample + phase);
    phase = 2*pi*rand;                                      % 随机相位(非相干慢跳频)                          
%     phase = phase + 2*pi*f2*(fsample*Tc)/fsample;          % 连续相位(相干慢跳频)
end
FH_sin(signal_len+1:end) = []; 
%%%%%%%%%%%%%%%%%%%   跳频信号产生   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
St_mix_fh = FH_sin.*St_baseband;                            % 跳频调制
N0 = 48;                                                    % 宽带带通滤波器阶数
tr_fir_coff = fir1(N0,[2*fh/fsample,2*(fh+2*fm)/fsample]);%[0.25,0.5]);                        % 发射宽带带通滤波器系数
%%%%%%%%%%%%%%%%%%%    发射信号通过宽带带通滤波器   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
St = zeros(signal_len,1);                    % 通过宽带滤波器后的发射信号
for i = 1:N0/2
    St(i) = tr_fir_coff * [zeros(N0/2-i+1,1);St_mix_fh(1:N0/2+i)];
end
for i = (N0/2+1):(signal_len - N0/2)
    St(i) = tr_fir_coff * St_mix_fh(i-N0/2:N0/2+i);
end
for i = (signal_len - N0/2 + 1):signal_len
    St(i) = tr_fir_coff * [St_mix_fh(i-N0/2:signal_len);zeros(i+N0/2-signal_len,1)];
end
St = (St-mean(St))/std(St);                                   % 发射信号幅度归一化

figure
subplot(311),plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(St_baseband))),grid on,title('DPSK基带信号频谱');
subplot(312),plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(St_mix_fh))),grid on,title('跳频混频后信号频谱');
subplot(313),plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(St))),grid on,title('发射信号频谱');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%        信道,噪声及干扰设置       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 信道为自由空间传播信道 , 噪声为加性高斯白噪声
SNR = 10;                                    % 加性高斯白噪声信噪比dB
SJR = -5;                                    % 信干比 : 接收信号功率/干扰信号功率

snr = SNR - 10*log10(fsample/2/Bt);           % 折算为带内白噪声强度
Jam = zeros(signal_len,1);                  % 干扰信号
%%%%%%%%%%%%%%%%%%%    干扰信号产生   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Jam_type = 4;               % 1:单音干扰
                            % 2:多音干扰
                            % 3:噪声调频干扰
                            % 4:MFSK 干扰
                            % 其它:仅高斯白噪声干扰
                            
%%% 网控信道干扰率设置  %%%%%%%%%  
Jam_percent = 1/1;          % 干扰网控信道数比例
                            % 1/1: 干扰全部网控信道
                            % 1/2: 干扰1/2网控信道
                            % 1/3: 干扰1/3网控信道

Jam_ch_num = ceil(K*Jam_percent);        % 干扰网控信道数 
channel = randperm(K);
channel(Jam_ch_num+1:end) = [];

multi_sin_num = 4;                      % 多音干扰每信道干扰正弦波数
Jam_M = 4;                              % MFSK 干扰进制数


switch Jam_type
    case 1          %%%%%%%%%%% 单音干扰 %%%%%%%%%%%%%%
        f_jam = fh + fm + (channel'-K/2-1/2) * Bs;
        Jam = (cos(2*pi*f_jam*(1:signal_len)/fsample+2*pi*rand))';
        Jam = mean(Jam,2);
        Jam = (Jam - mean(Jam))/std(Jam);
        figure,plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Jam))),grid on,title('单音干扰信号频谱');
    case 2          %%%%%%%%%%% 多音干扰 %%%%%%%%%%%%%%     
        Jam = zeros(signal_len,multi_sin_num);
        for i = 1:Jam_ch_num
            f_jam = fh + fm + (channel(i)-K/2) * Bs - rand(multi_sin_num,1) * Bs;
            Jam = Jam + (cos(2*pi*f_jam*(1:signal_len)/fsample+2*pi*rand))';            
        end
        Jam = mean(Jam,2);
        Jam = (Jam - mean(Jam))/std(Jam);
        figure,plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Jam))),grid on,title([num2str(multi_sin_num),'音干扰信号频谱']);
    case 3          %%%%%%%%%%% 噪声调频干扰 %%%%%%%%%%
        for i = 1:Jam_ch_num
            f_jam = rand(signal_len,1)*Bs - Bs/2;
            Jam = Jam + amod(f_jam,fh+fm+(channel(i)-K/2-1/2)*Bs,fsample,'fm');
        end
        Jam = (Jam - mean(Jam))/std(Jam);
        figure,plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Jam))),grid on,title('噪声调频干扰信号频谱');
    case 4          %%%%%%%%%%%  MFSK 干扰   %%%%%%%%%%    
        for i = 1:Jam_ch_num
            Jam = Jam + DMOD(randsrc(Data_len,1,(0:Jam_M-1)),fm+fh+(channel(i)-K/2-1/2)*Bs, fs, fsample,'fsk', Jam_M, Bs/Jam_M);
        end
        Jam = (Jam - mean(Jam))/std(Jam);
        figure,plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Jam))),grid on,title([num2str(Jam_M),'FSK干扰信号频谱']);
    otherwise       %%%%%%%%%%%  无人为干扰  %%%%%%%%%%
        
end
Jam = Jam*10^(-SJR/20);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%         跳频信号接收             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% 设接收系统与接收信号严格同步 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rt = awgn(St,snr) + Jam;                              % 接收信号
figure
plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Rt))),grid on,title('接收信号频谱');

%%%   接收参数设置
N1 = 48;                                            % 宽带带通滤波器阶数
wide_fir_coff = fir1(N1,[2*fh/fsample,2*(fh+2*fm)/fsample]);                % 接收机宽带带通滤波器系数
N2 = 48;                                              % 宽带低通滤波器阶数
lowpass_fir_coff = fir1(N2,4*fm/fsample);             % 解跳后宽带低通滤波器系数
N3 = 96;                                              % 窄带带通滤波器阶数
bandpass_fir_coff = fir1(N3,[0.22,0.28]);             % 降采样后窄带带通滤波器系数
%%%%%%%%%%%%%%%%%%%    接收信号通过宽带带通滤波器   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rt_passband = zeros(signal_len,1);                    % 通过宽带滤波器后的信号
for i = 1:N1/2
    Rt_passband(i) = wide_fir_coff * [zeros(N1/2-i+1,1);Rt(1:N1/2+i)];
end
for i = (N1/2+1):(signal_len - N1/2)
    Rt_passband(i) = wide_fir_coff * Rt(i-N1/2:N1/2+i);
end
for i = (signal_len - N1/2 + 1):signal_len
    Rt_passband(i) = wide_fir_coff * [Rt(i-N1/2:signal_len);zeros(i+N1/2-signal_len,1)];
end

% figure
% plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Rt_passband))),grid on,title('接收信号宽带滤波后频谱');

%%%%%%%%%%%%%%%%%%%    本地跳频载波产生   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % FH_sin = zeros(signal_len,1);                               % 本地跳频信号载波
% % phase = 2*pi*rand;
% % for i = 1:ceil(signal_len/fsample/Tc)
% %     f2 = fh + (FH_seq(i)*2 - (K-1))*Bs/2;                   % 跳频频点中心频率
% %     FH_sin((i-1)*fsample*Tc + 1:i*fsample*Tc) = cos(2*pi*f2*(0:fsample*Tc-1)/fsample + phase);
% %     phase = 2*pi*rand;                                      % 随机相位(非相干慢跳频)                          
% % %     phase = phase + 2*pi*f*(fsample*Tc)/fsample;          % 连续相位(相干慢跳频)
% % end
% % FH_sin(signal_len+1:end) = []; 
%%%%%%%%%%%%%%%%%%%    解跳   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rt_mix_fh = Rt_passband.*FH_sin;                      % 接收信号与本地跳频信号混频

figure
plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Rt_mix_fh))),grid on,title('解跳混频后信号频谱');

%%%%%%%%%%%%%%%%%%%    解跳信号通过宽带低通滤波器   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rt_lowpass = zeros(signal_len,1);                    % 通过宽带低通滤波器后的信号
for i = 1:N2/2
    Rt_lowpass(i) = lowpass_fir_coff * [zeros(N2/2-i+1,1);Rt_mix_fh(1:N2/2+i)];
end
for i = (N2/2+1):(signal_len - N2/2)
    Rt_lowpass(i) = lowpass_fir_coff * Rt_mix_fh(i-N2/2:N2/2+i);
end
for i = (signal_len - N2/2 + 1):signal_len
    Rt_lowpass(i) = lowpass_fir_coff * [Rt_mix_fh(i-N2/2:signal_len);zeros(i+N2/2-signal_len,1)];
end
figure
plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Rt_lowpass))),grid on,title('解跳混频去除倍频信号频谱');

%%%%%%%%%%%%%%%%%%%    采样   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sample_rate = fm*8;                                   % 8 倍载波采样速率
Rt_bassband = Rt_lowpass(1:fsample/sample_rate:end);  % 基带MFSK信号(降采样)
sample_len = length(Rt_bassband);
figure
plot((0:sample_len-1)/(sample_len)*sample_rate,abs(fft(Rt_bassband))),grid on,title('基带采样信号频谱');
%%%%%%%%%%%%%%%%%%%    窄带滤波   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rt_baseband = zeros(sample_len,1);
for i = 1:N3/2
    Rt_baseband(i) = bandpass_fir_coff * [zeros(N3/2-i+1,1);Rt_bassband(1:N3/2+i)];
end
for i = (N3/2+1):(sample_len - N3/2)
    Rt_baseband(i) = bandpass_fir_coff * Rt_bassband(i-N3/2:N3/2+i);
end
for i = (signal_len - N3/2 + 1):sample_len
    Rt_baseband(i) = bandpass_fir_coff * [Rt_bassband(i-N3/2:sample_len);zeros(i+N3/2-sample_len,1)];
end
figure
plot((0:sample_len-1)/(sample_len)*sample_rate,abs(fft(Rt_baseband))),grid on,title('窄带滤波后基带信号频谱');

%%%%%%%%%%%%%%%%%%    DPSK 信号解调    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rdata_diff = DDEMOD(Rt_baseband, fm, fs, sample_rate, 'psk/noncoherence',2);
Rdata = xor(Rdata_diff(1:end-1),Rdata_diff(2:end));                 % 接收数据
error = length(find(Data ~= Rdata))/length(Data)*100;               % 误码率
figure
subplot(311),plot(Data),grid on,title('发送数据');
subplot(312),plot(Rdata),grid on,title('接收数据');
subplot(313),plot(Data-Rdata),grid on,title(['误码率=',num2str(error),'%']);