clc,clear
close all

% 慢跳頻系統(tǒng)
% 數(shù)據(jù)采用 MFSK 調(diào)制


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%         跳頻參數(shù)設置             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = 4;              % MFSK
T = 0.001;          % 數(shù)據(jù)比特周期 (s)

L = log2(M);
Ts = L*T;            % MFSK 的碼符號周期
fm = 50000;          % MFSK 信號載頻中心頻率
fs = 1/Ts;           % MFSK 的 M 個單載頻的頻率間隔 (fs >= 1/Ts)
Bs = M*fs;           % MFSK 信號帶寬 (跳頻頻率最小間隔)

Tc = 4*Ts;           % 跳頻周期 (Tc>=Ts)
K = 32;              % 跳頻頻率個數(shù)
fh = 100000;         % 跳頻信號中心頻率
Bt = K*Bs;          % 發(fā)射信號帶寬

t = 0.1;              % 仿真時間
fsample = 800000;   % 采樣速率 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%         跳頻信號產(chǎn)生             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Data_len = floor(t/Ts);                                     % 用戶數(shù)據(jù)長度
Data = randsrc(Data_len,1,(0:M-1));                         % 用戶數(shù)據(jù)
%%%%%%%%%%%%%%%%        MFSK 調(diào)制       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

St_baseband = DMOD(Data,fm, fs, fsample,'fsk', M);         % MFSK 信號
% figure, plot(St_baseband),title('MFSK基帶信號');
signal_len = length(St_baseband);                           % 仿真時間內(nèi)采樣數(shù)據(jù)長度
%%%%%%%%%%%%%%%%%%%   跳頻載波產(chǎn)生   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;          % 連續(xù)相位(相干慢跳頻)
end
FH_sin(signal_len+1:end) = []; 
%%%%%%%%%%%%%%%%%%%   跳頻信號產(chǎn)生   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
St_mix_fh = FH_sin.*St_baseband;                            % 跳頻調(diào)制
N0 = 48;                                                    % 寬帶帶通濾波器階數(shù)
tr_fir_coff = fir1(N0,[2*fh/fsample,2*(fh+2*fm)/fsample]);%[0.25,0.5]);                        % 發(fā)射寬帶帶通濾波器系數(shù)
%%%%%%%%%%%%%%%%%%%    發(fā)射信號通過寬帶帶通濾波器   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
St = zeros(signal_len,1);                    % 通過寬帶濾波器后的發(fā)射信號
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);                                   % 發(fā)射信號幅度歸一化

figure
subplot(311),plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(St_baseband))),grid on,title([num2str(M),'FSK基帶信號頻譜']);
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('發(fā)射信號頻譜');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%        信道,噪聲及干擾設置       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 信道為自由空間傳播信道 , 噪聲為加性高斯白噪聲
SNR = 0;                                    % 加性高斯白噪聲信噪比dB
SJR = -5;                                    % 信干比 : 接收信號功率/干擾信號功率

snr = SNR - 10*log10(fsample/2/Bt);           % 折算為帶內(nèi)白噪聲強度
Jam = zeros(signal_len,1);                  % 干擾信號
%%%%%%%%%%%%%%%%%%%    干擾信號產(chǎn)生   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Jam_type = 4;               % 1:單音干擾
                            % 2:多音干擾
                            % 3:噪聲調(diào)頻干擾
                            % 4:MFSK 干擾
                            % 其它:僅高斯白噪聲干擾
                            
%%% 網(wǎng)控信道干擾率設置  %%%%%%%%%  
Jam_percent = 1/2;          % 干擾網(wǎng)控信道數(shù)比例
                            % 1/1: 干擾全部網(wǎng)控信道
                            % 1/2: 干擾1/2網(wǎng)控信道
                            % 1/3: 干擾1/3網(wǎng)控信道
                            
Jam_ch_num = ceil(K*Jam_percent);        % 干擾網(wǎng)控信道數(shù) 
channel = randperm(K);
channel(Jam_ch_num+1:end) = [];

multi_sin_num = 4;                      % 多音干擾每信道干擾正弦波數(shù)
Jam_M = 4;                              % MFSK 干擾進制數(shù)


switch Jam_type
    case 1          %%%%%%%%%%% 單音干擾 %%%%%%%%%%%%%%
        f_jam = fh + fm + (channel'-K/2-1/2) * Bs + randsrc(Jam_ch_num,1,(1-M/2:M/2)) * fs - fs/2;
        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          %%%%%%%%%%% 噪聲調(diào)頻干擾 %%%%%%%%%%
        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('噪聲調(diào)頻干擾信號頻譜');
    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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%         跳頻信號接收             %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%% 設接收系統(tǒng)與接收信號嚴格同步 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rt = awgn(St,snr) + Jam;                              % 接收信號
figure
plot((0:signal_len-1)/(signal_len)*fsample,abs(fft(Rt))),grid on,title('接收信號頻譜');

%%%   接收參數(shù)設置
N1 = 48;                                            % 寬帶帶通濾波器階數(shù)
wide_fir_coff = fir1(N1,[2*fh/fsample,2*(fh+2*fm)/fsample]);                % 接收機寬帶帶通濾波器系數(shù)
N2 = 48;                                              % 寬帶低通濾波器階數(shù)
lowpass_fir_coff = fir1(N2,4*fm/fsample);             % 解跳后寬帶低通濾波器系數(shù)
N3 = 96;                                              % 窄帶帶通濾波器階數(shù)
bandpass_fir_coff = fir1(N3,[0.22,0.28]);             % 降采樣后窄帶帶通濾波器系數(shù)
%%%%%%%%%%%%%%%%%%%    接收信號通過寬帶帶通濾波器   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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('接收信號寬帶濾波后頻譜');

%%%%%%%%%%%%%%%%%%%    本地跳頻載波產(chǎn)生   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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;          % 連續(xù)相位(相干慢跳頻)
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('窄帶濾波后基帶信號頻譜');

%%%%%%%%%%%%%%%%%%    MFSK 信號解調(diào)    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rdata = DDEMOD(Rt_baseband, fm, fs, sample_rate, 'fsk/noncoherence',M);
error = length(find(Data ~= Rdata))/length(Data)*100;           % 誤碼率
figure
subplot(311),plot(Data),grid on,title('發(fā)送數(shù)據(jù)');
subplot(312),plot(Rdata),grid on,title('接收數(shù)據(jù)');
subplot(313),plot(Data-Rdata),grid on,title(['誤碼率 = ',num2str(error),'%']);