通信干扰识别_single sinusoid jam

一、目标

• 生成单音干扰、多音干扰、宽带噪声干扰、部分频带噪声干扰、宽带梳状谱干扰、线性调频干扰等6 种通信干扰信号。
• 选择合适的特征参数，采用决策树法实现对上述干扰信号的识别，高斯白噪声信道，干噪比(JNR)为0~15dB，识别正确率大于95%。
• 选择合适的特征参数，采用NN 或者SVM 机器学习实现对上述干扰信号的识别，高斯白噪声信道，干噪比(JNR)为0~15dB，识别正确率大于95%。

二、通信干扰信号生成

(一) 单音干扰

单音干扰信号在某个频点上发射，是一个单频连续正弦波形。时域表达式如下
$$J(t)=\sqrt{P_J}\cos (2\pi f_{J}t+{\theta }_J),{\theta }_J\sim (0,2\pi )$$
其中$P_J$为干扰信号功率，$f_J$为干扰信号频率，${\theta }_J$为干扰信号相位且服从均匀分布。

% produce a fragment of single sinusoid jam
clear;close all;
Pl = 1;
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
fs = 50;
x = sin(2*pi*fs*t); 
jnr = 200;
y = awgn(x,jnr,'measured');   % Sinusoids plus noise   
figure;set(gcf,'Color','w','Position',[400 300 600 300]);
plot(Fs*t(1:100),y(1:100),'Color',[0,0.4,0.8]);grid on;
title('single-tone jam')
xlabel('time (milliseconds)');
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure;set(gcf,'Color','w','Position',[400 300 600 300]);
semilogy(f,2*abs(Y(1:NFFT/2+1)),'Color',[0,0.4,0.8]);grid on;
title('Single-Sided Amplitude Spectrum of single-tone jam');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
ylim([10e-4,1]);grid on;
% C par
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

(二) 多音干扰

多音干扰信号在多个离散的频点上发射，时域表达式如下
$$J(t)=\sum _{l=1}^Q\sqrt{P_l}\cos (2\pi f_{l}t+{\theta }_l),{\theta }_l\sim U(0,2\pi )$$


由Figure 4可以看出，多音干扰信号从100Hz的频点开始，每间隔10Hz有一个冲激，代表一个干扰频点，一个有10个冲激。

% produce a fragment of multiple sinusoid jam
clear;close all;jnr = 200; % jam-noise ratio
Q = 10; % number of multiple sinusoid
Pl = 1-0.1*rand(1,Q);                       % Jam power

Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid
fs0 = 100;deltafs = 10;
fs = fs0:deltafs:fs0+(Q-1)*deltafs; % produce same-distance fs
theta = 2*pi*rand(1,Q);% random phase
xx = cos(2*pi*fs'*t+theta'*ones(1,length(t)));% cos
x = sqrt(Pl)*xx; % multi frequency Sinusoids jam with random phase

x_mid = (max(x)+min(x))/2;
x_nor = 2*(x-x_mid)/(max(x)-min(x));% normalize the data into [-1,1]
x_dec = (x_nor-mean(x_nor))/sqrt(var(x_nor));% decentralize the data

figure;
subplot(2,1,1);
plot(Fs*t(1:5*Fs/fs(1)),x(1:5*Fs/fs(1)))
title('original data')
xlabel('time (milliseconds)');grid on;
subplot(2,1,2);
plot(Fs*t(1:5*Fs/fs(1)),x_dec(1:5*Fs/fs(1)))
title('pre_processed data')
xlabel('time (milliseconds)');grid on;

y = awgn(x,jnr,'measured');    % Sinusoids plus noise

figure;set(gcf,'Color','w','Position',[400 300 600 300]);
plot(Fs*t(1:10*Fs/fs(1)),y(1:10*Fs/fs(1)),'Color',[0,0.4,0.8]);grid on;
title('multiple-tone jam')
xlabel('time (milliseconds)');
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
figure;set(gcf,'Color','w','Position',[400 300 600 300]);
semilogy(f,2*abs(Y(1:NFFT/2+1)),'Color',[0,0.4,0.8]);grid on;
title('Single-Sided Amplitude Spectrum of multiple-tone jam');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
ylim([10e-4,1]);grid on;

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48

(三) 宽带噪声干扰

宽带噪声干扰是将高斯白噪声通过一个宽带滤波器，得到在每一个宽带范围内的阻塞噪声。设白噪声$\epsilon (t)\tilde{N}(0,P_J/2W_I)$ ，其中$P_J$是干扰功率，$W_I$是低通滤波器的截止频率。则宽带噪声干扰可以表示为 $n(t)=\epsilon (t)\otimes h(t)$，滤波器的频域表达式为
H ( j 2 π f ) = { 1 , ∣ f ∣ ≤ W I 0 , o t h e r w i s e H(j2\pi f) = \left\{

$$\begin{array}{l} 1,{\rm{ }}|f| \le {W_I}\\ 0,{\rm{ }}otherwise \end{array}$$ \right. H(j2πf)={1,∣f∣≤WI​0,otherwise​

% produce a fragment of narrow band noise 
clear;close all;jnr = 20; % jam-noise ratio
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal
t = (0:L-1)*T;                % Time vector
P=1;% power of noise
y = wgn(L,1,P); % white noise
% 
% y = awgn(x,jnr,'measured');    % Sinusoids plus noise
figure;plot(1000*t(1:L/2),y(1:L/2));
title('White Noise')
xlabel('time (milliseconds)');grid on;

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure;plot(f,2*10*log10(abs(Y(1:NFFT/2+1)))) 
title('Single-Sided Amplitude Spectrum of White Noise')
xlabel('Frequency (Hz)');ylabel('|Y(f)|(dB)');grid on;
% ylim([10e-4,1])

%% bandpass filter
WI = 50; % window of narrow band
FJ = 600;% jamming frequency
fl_kaiser = [FJ-WI/2 FJ-WI/4 FJ+WI/4 FJ+WI/2];
fl_mag = [0 1 0];
fl_dev = [0.05 0.05 0.05];
[fl_n_kaiser,fl_wn,fl_beta,fl_ftype]=kaiserord(fl_kaiser,fl_mag,fl_dev,Fs);
h= fir1(fl_n_kaiser,fl_wn,fl_ftype,kaiser(fl_n_kaiser+1,fl_beta));
Y_bp = filter(h,1,y);% the result of filter in time domain
figure;
freqz(h)
figure;
subplot(2,1,1);plot(1000*t(1:L/2),y(1:L/2));% original noise
title('White Noise');xlabel('time (milliseconds)');grid on;
subplot(2,1,2);plot(1000*t(1:L/2),Y_bp(1:L/2));% original noise
title('Narrow Band White Noise');xlabel('time (milliseconds)');grid on;

figure;set(gcf,'Color','w','Position',[400 300 600 300]);
plot(1000*t(1:L/2),Y_bp(1:L/2),'Color',[0,0.4,0.8]);grid on;
title('Narrowband Noise jam')
xlabel('time (milliseconds)');


NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(Y_bp,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure;set(gcf,'Color','w','Position',[400 300 600 300]);
plot(f,2*10*log10(abs(Y(1:NFFT/2+1))),'Color',[0,0.4,0.8]);grid on;
title('Single-Sided Amplitude Spectrum of Narrowband Noise jam')
xlabel('Frequency (Hz)');ylabel('|Y(f)|(dB)');grid on;
% ylim([10e-4,1])


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57

(四) 窄带噪声干扰

与宽带噪声干扰产生

方法类似，将白噪声通过一个窄带滤波器，得到瞄准某一个窄带范围的阻塞信号。窄带滤波器的频域表达式为
H ( j 2 π f ) = { 1 , ∣ f − f J ∣ ≤ W I 2 0 , o t h e r w i s e H(j2\pi f) = \left\{

$$\begin{array}{l} 1,{\rm{ }}|f - {f_J}| \le \frac{ { {W_I}}}{2}\\ 0,{\rm{ }}otherwise \end{array}$$ \right. H(j2πf)={1,∣f−fJ​∣≤2WI​​0,otherwise​
$f_J$表示瞄准频率。设定$f_J=600\text{ Hz}$ ，$W_I=50\text{ Hz}$ ，得到时域窄带噪声干扰信号和其单边频谱如图Figure 7、Figure 8所示。

% produce a fragment of broad band noise 
clear;close all;
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal
t = (0:L-1)*T;                % Time vector
P=1;% power of noise
y = wgn(L,1,P/600); % white noise
% 
% y = awgn(x,jnr,'measured');    % Sinusoids plus noise
figure;plot(1000*t(1:L/2),y(1:L/2));
title('White Noise')
xlabel('time (milliseconds)');grid on;

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure;plot(f,2*10*log10(abs(Y(1:NFFT/2+1)))) 
title('Single-Sided Amplitude Spectrum of White Noise')
xlabel('Frequency (Hz)');ylabel('|Y(f)|(dB)');grid on;
% ylim([10e-4,1])

%% bandpass filter
WI = 600; % window of broad band
fl_kaiser = [WI WI+50];
fl_mag = [1 0];
fl_dev = [0.05 0.01];
[fl_n_kaiser,fl_wn,fl_beta,fl_ftype]=kaiserord(fl_kaiser,fl_mag,fl_dev,Fs);
h= fir1(fl_n_kaiser,fl_wn,fl_ftype,kaiser(fl_n_kaiser+1,fl_beta));
Y_bp = filter(h,1,y);% the result of filter in time domain
figure;freqz(h);
figure;
subplot(2,1,1);plot(1000*t(1:L/2),y(1:L/2));% original noise
title('White Noise');xlabel('time (milliseconds)');grid on;
subplot(2,1,2);plot(1000*t(1:L/2),Y_bp(1:L/2));% original noise
title('Broad Band White Noise');xlabel('time (milliseconds)');grid on;

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(Y_bp,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure;plot(f,2*10*log10(abs(Y(1:NFFT/2+1)))) 
title('Single-Sided Amplitude Spectrum of Narrow band White Noise')
xlabel('Frequency (Hz)');ylabel('|Y(f)|(dB)');grid on;
% ylim([10e-4,1])



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

(五) 宽带梳状谱干扰

将白噪声通过多个窄带滤波器，得到瞄准多个窄带范围的阻塞信号。宽带梳状谱干扰的滤波器频域表达式为
H ( j 2 π f ) = { 1 , ∣ f − f J , i ∣ ≤ W I , i 2 0 , o t h e r w i s e H(j2\pi f) = \left\{

$$\begin{array}{l} 1,{\rm{ }}|f - {f_{J,i}}| \le \frac{ { {W_{I,i}}}}{2}\\ 0,{\rm{ }}otherwise \end{array}$$ \right. H(j2πf)={1,∣f−fJ,i​∣≤2WI,i​​0,otherwise​
其中 $i=1,2,\dots ,Q$，$Q$ 表示窄带滤波器的个数。设定 $Q=4$，第$i$个窄带滤波器的瞄准频率$f_{J,i}=100\mathrm{Hz}+(i-1)\times 100\mathrm{Hz},i=1,2,...,Q$ ， $W_{I,i}=50\mathrm{Hz}$。得到时域宽带梳状谱干扰信号和其单边频谱如图Figure 9、Figure 10所示。

% produce a fragment of comb-like noise 
clear;close all;
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal
t = (0:L-1)*T;                % Time vector
P=1;% power of noise
y = wgn(L,1,P); % white noise
% 
% y = awgn(x,jnr,'measured');    % Sinusoids plus noise
figure;plot(1000*t(1:L/2),y(1:L/2));
title('White Noise')
xlabel('time (milliseconds)');grid on;

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure;plot(f,2*10*log10(abs(Y(1:NFFT/2+1)))) 
title('Single-Sided Amplitude Spectrum of White Noise')
xlabel('Frequency (Hz)');ylabel('|Y(f)|(dB)');grid on;
% ylim([10e-4,1])

%% bandpass filter
WI = 50; % window of broad band
FJ_set = [100,200,300,400];
for i = 1:4
    FJ = FJ_set(i);% jamming frequency
    fcuts = [FJ-WI/3 FJ-WI/4 FJ+WI/4 FJ+WI/3];
    mags = [0 1 0];
    devs = [0.02 0.02 0.02];
    [n,wn,beta,ftype]=kaiserord(fcuts,mags,devs,Fs);
    h= fir1(n,wn,ftype,kaiser(n+1,beta));
    Y_bp_temp(i,:) = filter(h,1,y);% the result of filter in time domain
end
Y_bp = sum(Y_bp_temp);
figure;
subplot(2,1,1);plot(1000*t(1:L/2),y(1:L/2));% original noise
title('White Noise');xlabel('time (milliseconds)');grid on;
subplot(2,1,2);plot(1000*t(1:L/2),Y_bp(1:L/2));% original noise
title('Partialband Noise Jam');xlabel('time (milliseconds)');grid on;

figure;set(gcf,'Color','w','Position',[400 300 600 300]);
plot(1000*t(1:L/2),Y_bp(1:L/2),'Color',[0,0.4,0.8]);grid on;
title('Partialband Noise Jam')
xlabel('time (milliseconds)');

NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(Y_bp,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure;set(gcf,'Color','w','Position',[400 300 600 300]);
plot(f,2*10*log10(abs(Y(1:NFFT/2+1))),'Color',[0,0.4,0.8]);grid on;
title('Single-Sided Amplitude Spectrum of Partialband Noise Jam')
xlabel('Frequency (Hz)');ylabel('|Y(f)|(dB)');grid on;
% ylim([10e-4,1])


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

(六) 扫频干扰

扫频干扰在部分频谱内呈现出随时间线性扫描的特性。其时域表达式为
$$J(t)=A(t)\exp (j{\beta }_i{t}^2/2+j{\omega }_{i}t+j\phi )$$
A ( t ) = A r e c t ( t / T ) { A , ∣ t ∣ ≤ T / 2 0 , e l s e A(t) = Arect(t/T)\left\{

$$\begin{array}{l} A,{\rm{ }}|t| \le T/2\\ 0,{\rm{ }}else \end{array}$$ \right. A(t)=Arect(t/T){A,∣t∣≤T/20,else​
其中${\beta }_i$为扫频速率，${\omega }_i$为初始角频率， $\phi$为初始相位，$T$为扫频持续时间。由时域表达式可以看出，扫频干扰的瞬时频率随时间呈不断的线性变化。设定${\beta }_i=100\pi rad/s$ ，${\omega }_i=20\pi \mathrm{rad}$ ，$\phi =0$ ，$T=2\mathrm{s}$ ，得到扫频干扰时域信号和其单边频谱如图Figure 11、 Figure 12所示。

% produce a fragment of scanning frequency jam
clear;close all;jnr = 1000; % jam-noise ratio
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal 2Second
t = (0:L-1)*T;                % Time vector
om = 2*pi*10;
be = 2*pi*50;
phi = 0;
x = exp(1i*0.5*be*t.^2+1i*om*t+1i*phi);
% x = exp(1i*om*t+1i*phi);

x = real(x);
y = awgn(x,jnr,'measured');    % Sinusoids plus noise
figure;plot(Fs*t,y)
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)');grid on;

figure;set(gcf,'Color','w','Position',[400 300 600 300]);
plot(1000*t(1:L/2),y(1:L/2),'Color',[0,0.4,0.8]);grid on;
title('Sweeping Jam')
xlabel('time (milliseconds)');


NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
figure;set(gcf,'Color','w','Position',[400 300 600 300]);
plot(f,2*10*log10(abs(Y(1:NFFT/2+1))),'Color',[0,0.4,0.8]);grid on;
title('Single-Sided Amplitude Spectrum of Sweeping Jam')
xlabel('Frequency (Hz)');ylabel('|Y(f)|(dB)');grid on;
xlim([0,500])

% plot(2*abs(Y(1:NFFT/2+1)))



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

干扰信号产生程序

% produce the jamming signal
function x = f_produce_jam(jnr,typein,type)
% type = 1 静态参数生成 type=2 动态生成
switch typein
    case 1
        x = f_single_jam(jnr,type); % test signal
    case 2
        x = f_multi_jam(jnr,type); % test signal
    case 3
        x = f_narrowbandnoise_jam(jnr,type); % test signal
    case 4
        x = f_broadbandnoise_jam(jnr,type); % test signal
    case 5
        x = f_combnoise_jam(jnr,type); % test signal
    case 6
        x = f_scanning_jam(jnr,type); % test signal
end

% produce a fragment of single sinusoid jam
function y = f_single_jam(jnr,type)
if type == 1   % 静态参数
    fs = 50;
elseif  type == 2     % 动态参数
    fs = randi([50,600]);
end
Pl = 1;
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

x = Pl*sin(2*pi*fs*t); 
y = awgn(x,jnr,'measured');   % Sinusoids plus noise

% produce a fragment of multiple sinusoid jam
function y=f_multi_jam(jnr,type)
if type == 1   % 静态参数
    Q = 10; % number of multiple sinusoid
    fs0 = 100; % start freq
elseif  type == 2     % 动态参数
    Q = randi([10,15]); % number of multiple sinusoid dynamic produce
    fs0 = randi([50,350]); 
end

Pl = 1-0.1*rand(1,Q);                       % Jam power
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid

deltafs = 40;
fs = fs0:deltafs:fs0+(Q-1)*deltafs; % produce same-distance fs
theta = 2*pi*rand(1,Q);% random phase
xx = cos(2*pi*fs'*t+theta'*ones(1,length(t)));% cos
x = sqrt(Pl)*xx; % multi frequency Sinusoids jam with random phase
y = awgn(x,jnr,'measured');    % Sinusoids plus noise

% produce a fragment of broad band noise 
function y=f_broadbandnoise_jam(jnr,type)
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal
P=1;% power of noise
y = wgn(L,1,P); % white noise

%% bandpass filter
if type == 1   % 静态参数
    WI = 800; % window of broad band
elseif  type == 2     % 动态参数
    WI = randi([800,900]); 
end

fl_kaiser = [WI WI+50];
fl_mag = [1 0];
fl_dev = [0.05 0.01];
[fl_n_kaiser,fl_wn,fl_beta,fl_ftype]=kaiserord(fl_kaiser,fl_mag,fl_dev,Fs);
h= fir1(fl_n_kaiser,fl_wn,fl_ftype,kaiser(fl_n_kaiser+1,fl_beta));
Y_bp = filter(h,1,y);% the result of filter in time domain
y = awgn(Y_bp,jnr,'measured');    %  plus noise


% produce a fragment of narrow band noise
function y=f_narrowbandnoise_jam(jnr,type)

Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal
P=1;% power of noise
y = wgn(L,1,P); % white noise
%% bandpass filter
WI = 50; % window of narrow band

if type == 1   % 静态参数
    FJ = 400;% jamming frequency
elseif  type == 2     % 动态参数
    FJ = randi([50,600]); 
end
fl_kaiser = [FJ-WI/2 FJ-WI/4 FJ+WI/4 FJ+WI/2];
fl_mag = [0 1 0];
fl_dev = [0.05 0.05 0.05];
[fl_n_kaiser,fl_wn,fl_beta,fl_ftype]=kaiserord(fl_kaiser,fl_mag,fl_dev,Fs);
h= fir1(fl_n_kaiser,fl_wn,fl_ftype,kaiser(fl_n_kaiser+1,fl_beta));
Y_bp = filter(h,1,y);% the result of filter in time domain
y = awgn(Y_bp,jnr,'measured');    %  plus noise

% produce a fragment of comb-like noise 
function y=f_combnoise_jam(jnr,type)
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal
P=1;% power of noise
y = wgn(L,1,P); % white noise

%% bandpass filter
if type == 1   % 静态参数
    WI = 50; % window of broad band
    FJ_set = [100,200,300,400];
elseif  type == 2     % 动态参数
    WI = randi([40,60]); % number of multiple sinusoid dynamic produce
    num = randi([4,5]);
    startf = randi([90,500]);
    freq_delta = randi([50,150]);
    FJ_set = startf:100:(num-1)*100+startf;
end

for i = 1:size(FJ_set,2)
    FJ = FJ_set(i);% jamming frequency
    fcuts = [FJ-WI/3 FJ-WI/4 FJ+WI/4 FJ+WI/3];
    mags = [0 1 0];
    devs = [0.02 0.02 0.02];
    [n,wn,beta,ftype]=kaiserord(fcuts,mags,devs,Fs);
    h= fir1(n,wn,ftype,kaiser(n+1,beta));
    Y_bp_temp(i,:) = filter(h,1,y);% the result of filter in time domain
end
Y_bp = sum(Y_bp_temp);
y = awgn(Y_bp,jnr,'measured');    %  plus noise
% produce a fragment of single sinusoid jam
function y = f_scanning_jam(jnr,type)
Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2*Fs;                     % Length of signal 2Second
t = (0:L-1)*T;                % Time vector
if type == 1   % 静态参数
    om = 2*pi*200;
    be = 2*pi*50;
    phi = 0;
elseif  type == 2     % 动态参数
    om = 2*pi*randi([0,500]); % number of multiple sinusoid dynamic produce
    be = 2*pi*randi([50,100]); 
    phi =2*pi*rand(1); 
end

x = exp(1i*0.5*be*t.^2+1i*om*t+1i*phi);
x = real(x);
y = awgn(x,jnr,'measured');   % Sinusoids plus noise


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

三、 特征参数提取

(一) 信号预处理

在对信号进行特征参数提取前，需要首先对信号进行预处理，包括归一化、中心化。

1. 归一化

设样本数据为$x(n),n=1,2,...,N$ ，采用中间值法将样本数据归一化到[-1,1]内，公式表示为
$$y(k)=\frac{x(k)-x_{\mathrm{mid}}}{\frac{1}{2}(x_{\max }-x_{\min })},k=1,2,...,N$$
$$x_{\mathrm{mid}}=\frac{x_{\max }+x_{\min }}{2}$$
其中 $x_{\max }$表示最大值，$x_{\min }$表示最小值。

2. 中心化

采用Z分法
$$y(k)=\frac{x(k)-\mu }{\sigma },k=1,2,...,N$$

(二) 时域参数

1. R参数

$$R=\frac{{\sigma }^2}{{\mu }^2}$$
其中$\mu$为随机变量$X$的均值，$\sigma $ 为其标准差，反映了变量$X$包络的变化程度。$R$越大，说明信号包络的变化程度越大；当$R$较小时，说明信号包络的变化程度较小。对6种类型的通信干扰信号计算其R参数，得到如Figure 13所示结果。

不同干扰信号R参数随JNR变化趋势
可以看出，单音干扰和扫频干扰由于其信号包络恒定，因此其R参数相比其他4种干扰信号要小，并且二者的R参数值几乎相同，随着JNR增大，R参数减小。宽带噪声干扰和窄带噪声干扰其信号包络与JNR大小无关，因此其R参数几乎不随JNR变化。梳状谱干扰随着JNR增大，R参数增大，说明其JNR越大时梳状谱信号的包络变化程度越大。多音干扰的R参数没有明显的特征。

% calculate R of different JNR
% plot the trend
clear;close all;
% 考虑包络 均值如何定义？
JNR = -5:1:20;
mm = 100;% Mont-carlo times 求平均？
type = 2;
for k=1:mm
    for i = 1:length(JNR)
        jnr = JNR(i);
        mj = f_multi_jam(jnr,type); % multilple freq
        R(1,i,k) = f_par_R(mj);
        sj = f_single_jam(jnr,type);% single freq
        R(2,i,k) = f_par_R(sj);
        nj = f_narrowbandnoise_jam(jnr,type); % narrow band noise jam
        R(3,i,k) = f_par_R(nj);
        bj = f_broadbandnoise_jam(jnr,type);% broad band noise jam
        R(4,i,k) = f_par_R(bj);
        cj = f_combnoise_jam(jnr,type);% comb-like noise jam
        R(5,i,k) = f_par_R(cj);  
        sj = f_scanning_jam(jnr,type);% comb-like noise jam
        R(6,i,k) = f_par_R(sj);  
    end
end
for  i = 1:6
    for j = 1:length(JNR)
        RR(i,j) = mean(R(i,j,:));
    end
end


row = size(RR,1);
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:row
    plot(JNR,RR(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('多音干扰','单音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('R参数')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

2. 时域矩偏度系数a3

$$a_3=\frac{E{(x-\mu )}^3}{{\sigma }^3}$$
$\mu$为时域信号包络 的均值，$\sigma$为其标准差。 $a_3$描述了频谱偏离正态分布的程度。$a_3$越大，说明频谱偏离正态分布的程度越大。

% calculate a3 of different JNR
% plot the trend
clear;close all;
JNR = -5:1:20;
mm = 1;% Mont-carlo times 求平均？
type = 2;
for k=1:mm
    for i = 1:length(JNR)
        jnr = JNR(i);
        mj = f_multi_jam(jnr,type); % multilple freq
        a3(1,i,k) = f_par_a3(mj);
        sj = f_single_jam(jnr,type);% single freq
        a3(2,i,k) = f_par_a3(sj);
        nj = f_narrowbandnoise_jam(jnr,type); % narrow band noise jam
        a3(3,i,k) = f_par_a3(nj);
        bj = f_broadbandnoise_jam(jnr,type);% broad band noise jam
        a3(4,i,k) = f_par_a3(bj);
        cj = f_combnoise_jam(jnr,type);% comb-like noise jam
        a3(5,i,k) = f_par_a3(cj);
        sj = f_scanning_jam(jnr,type);% comb-like noise jam
        a3(6,i,k) = f_par_a3(sj);     
    end
end
for  i = 1:5
    for j = 1:length(JNR)
        a3a(i,j) = mean(a3(i,j,:));
    end
end


row = size(a3a,1);
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:row
    plot(JNR,a3a(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('多音干扰','单音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('a3参数')
% Chracater R reflect the changing trend of waveform
function a3 = f_par_a3(x)
% Normalize the original data
x_mid = (max(x)+min(x))/2;
x_nor = 2*(x-x_mid)/(max(x)-min(x));% normalize the data into [-1,1]
x_dec = (x_nor-mean(x_nor))/sqrt(var(x_nor));% decentralize the data

x_en= x_dec.^2;%信号包络
sigma_x = sqrt(var(x_en));
mu_x = mean(x_en);
a3 = mean( (x_en-mu_x).^3)/(sigma_x^3);
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

(三) 频域参数

1. 载波因子系数C

对时域信号$x(n),n=1,2,...,N$进行FFT得到频域离散信号 $X[n],n=1,2,...,NFFT$，按照幅值大小进行排序得到$X[{\lambda }_1],X[{\lambda }_2],...,X[{\lambda }_{NFFT}]$ ，其信号谱最大值$X[{\lambda }_1]$与相邻较大值X[{\lambda _2}]之比，定义为载波因子系数$C$，描述了信号谱线的突出程度。
$$C=\frac{X[{\lambda }_1]}{X[{\lambda }_2]}$$

不同干扰信号C参数随JNR变化趋势

由于单音干扰其频域存在单一的一个类似冲激的信号，因此其载波因子系数C值较大，对于其他干扰信号，至少存在两个以上比较高的谱峰，因此载波因子系数C普遍较低。根据这一特征，可以将单音干扰信号与其他信号进行有效区别。

% calculate b3 of different JNR
% plot the trend
clear;close all;
JNR = -5:1:20;
mm = 1;% Mont-carlo times 求平均？
type = 2;% 动态
for k=1:mm
    for i = 1:length(JNR)
        jnr = JNR(i);
        mj = f_multi_jam(jnr,type); % multilple freq
        Y = f_cal_fft(mj);
        C(1,i,k) = f_par_C(Y);
        sj = f_single_jam(jnr,type);% single freq
        Y = f_cal_fft(sj);
        C(2,i,k) = f_par_C(Y);
        nj = f_narrowbandnoise_jam(jnr,type); % narrow band noise jam
        Y = f_cal_fft(nj);
        C(3,i,k) = f_par_C(Y);        
        bj = f_broadbandnoise_jam(jnr,type);% broad band noise jam
        Y = f_cal_fft(bj);
        C(4,i,k) = f_par_C(Y);        
        cj = f_combnoise_jam(jnr,type);% comb-like noise jam
        Y = f_cal_fft(cj);
        C(5,i,k) = f_par_C(Y);
        sj = f_scanning_jam(jnr,type);% comb-like noise jam
        Y = f_cal_fft(sj);
        C(6,i,k) = f_par_C(Y);   
    end
end
for  i = 1:6
    for j = 1:length(JNR)
        ca(i,j) = mean(C(i,j,:));
    end
end


row = size(ca,1);
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:row
    plot(JNR,ca(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('多音干扰','单音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('C 载波因子系数')
ylim([0,10])

% Chracater C reflect the changing trend of waveform
function C = f_par_C(Y)
%INPUT: single_wise spectrum
Y_s = sort(Y,'descend');
C = Y_s(1)/Y_s(20);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

2. 平均频谱平坦系数Fse

$$Fse=\sqrt{\frac{1}{N}\sum _{n=0}^{N_s-1}{(X_1(n)-\overline{X_1(n)})}^2}$$
$$X_1(n)=X(n)-X_2(n)$$
$$X_2(n)=\frac{1}{2L+1}\sum _{i=-L}^{L}X(n+i)$$
$X(n)$为频谱幅度，$\overline{X_1(n)}$为 均值，$X_1(n)$是对频谱进行平滑滤波。Fse反映了信号局部是否存在明显的冲激信号。

可以看出，单音干扰的Fse最大，多音次之，梳状谱、窄带、宽带依次减小，扫频干扰的Fse最低。Fse越大，意味着信号冲击部分波动较大，平坦度低，平坦系数较小时，冲激部分波动较小，平坦度高。

% calculate fse of different JNR
% plot the trend
clear;close all;
JNR = -5:1:20;
mm = 1;% Mont-carlo times 求平均？
type = 2;
for k=1:mm
    for i = 1:length(JNR)
        jnr = JNR(i);
        mj = f_multi_jam(jnr,type); % multilple freq
        Y = f_cal_fft(mj);
        Fse(1,i,k) = f_par_Fse(Y);
        sj = f_single_jam(jnr,type);% single freq
        Y = f_cal_fft(sj);
        Fse(2,i,k) = f_par_Fse(Y);
        nj = f_narrowbandnoise_jam(jnr,type); % narrow band noise jam
        Y = f_cal_fft(nj);
        Fse(3,i,k) = f_par_Fse(Y);        
        bj = f_broadbandnoise_jam(jnr,type);% broad band noise jam
        Y = f_cal_fft(bj);
        Fse(4,i,k) = f_par_Fse(Y);        
        cj = f_combnoise_jam(jnr,type);% comb-like noise jam
        Y = f_cal_fft(cj);
        Fse(5,i,k) = f_par_Fse(Y);     
        sj = f_scanning_jam(jnr,type);% comb-like noise jam
        Y = f_cal_fft(sj);
        Fse(6,i,k) = f_par_Fse(Y);   
    end
end
for  i = 1:6
    for j = 1:length(JNR)
        Fsea(i,j) = mean(Fse(i,j,:));
    end
end

row = size(Fsea,1);
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:2:3
    plot(JNR,Fsea(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('多音干扰','窄带干扰')
xlabel('JNR(dB)')
title('Fse 平均频谱平坦系数')
ylim([0,0.06])

% Chracater Fse reflect the changing trend of waveform
function Fse = f_par_Fse(Y)
% Y = Y/max(Y);% normalize the spectrum
P = Y;
L = floor(length(Y)*0.03);
for i = L+1:length(P)-L
    P1(i-L) = P(i)-sum(P(i-L:i+L))/(2*L+1);
end
Fse = sqrt(sum((P1-mean(P1)).^2)/length(P));
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

3. Rf参数

类似时域$R$参数的定义，给出频域 $R_f$参数的定义
$$R_f=\frac{{\sigma }^2}{{\mu }^2}$$
$\mu$为频谱幅度X的均值，$\sigma$为其标准差，反映了$X$包络的变化程度。$R_f$越大，说明信号的频谱包络变化程度越大，反之说明信号的频谱包络变化程度越小。

对不同干扰信号在不同JNR下计算其 参数，结果如图所示。可以看出，单音干扰、多音干扰以及窄带干扰的频域 参数值随JNR增大而逐渐增大，且单音干扰增大趋势较快。宽带干扰的$R_f$参数最小，扫频干扰的$R_f$参数也较低。在频谱上，宽带干扰和扫频干扰的幅度变化较小，也即其包络变化程度不够剧烈。据此，可以将多音干扰和窄带干扰与宽带干扰、梳状谱干扰以及扫频干扰进行区分。

% calculate Rf of different JNR
% plot the trend
clear;close all;
JNR = -5:1:20;
mm = 1;% Mont-carlo times 求平均？
type = 2;%  1 静态参数 2 动态参数 
for k=1:mm
    for i = 1:length(JNR)
        jnr = JNR(i);
        mj = f_multi_jam(jnr,type); % multilple freq
        Y = f_cal_fft(mj);
        Rf(1,i,k) = f_par_R_f(Y);
        sj = f_single_jam(jnr,type);% single freq
        Y = f_cal_fft(sj);
        Rf(2,i,k) = f_par_R_f(Y);
        nj = f_narrowbandnoise_jam(jnr,type); % narrow band noise jam
        Y = f_cal_fft(nj);
        Rf(3,i,k) = f_par_R_f(Y);        
        bj = f_broadbandnoise_jam(jnr,type);% broad band noise jam
        Y = f_cal_fft(bj);
        Rf(4,i,k) = f_par_R_f(Y);        
        cj = f_combnoise_jam(jnr,type);% comb-like noise jam
        Y = f_cal_fft(cj);
        Rf(5,i,k) = f_par_R_f(Y);        
        sj = f_scanning_jam(jnr,type);% comb-like noise jam
        Y = f_cal_fft(sj);
        Rf(6,i,k) = f_par_R_f(Y);   
    end
end
for  i = 1:6
    for j = 1:length(JNR)
        Rfa(i,j) = mean(Rf(i,j,:));
    end
end


row = size(Rfa,1);
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:row
    plot(JNR,Rfa(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('多音干扰','单音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('Rf 频域信号包络变化程度')
ylim([0,5])


% Chracater R in freq domain reflect the changing trend of waveform
function Rf = f_par_R_f(Y)
Y = Y/max(Y);
en = Y.^2;
uf = mean(en);
simgaf = sqrt(var(en));
Rf = simgaf.^2/uf.^2;
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

4. 频域矩偏度系数b3

$$b_3=\frac{E{(X-\mu )}^3}{{\sigma }^3}$$
$b_3$ 描述了频谱偏离正态分布的程度。$b_3$越大，说明频谱偏离正态分布的程度越大。

可以看出单音干扰偏离正态分布程度最大，其次是多音干扰，窄带干扰和梳状谱干扰。宽带干扰与正态分布类似，因为宽带干扰只是截断频带的高斯噪声，因此保留了近似随机的特性。

% calculate b3 of different JNR
% plot the trend
% clear;close all;
JNR = -5:1:20;
mm = 1;% Mont-carlo times 求平均？
type = 2;%  1 静态参数 2 动态参数 
for k=1:mm
    for i = 1:length(JNR)
        jnr = JNR(i);
        mj = f_multi_jam(jnr,type); % multilple freq
        Y = f_cal_fft(mj);
        b3(1,i,k) = f_par_b3(Y);
        sj = f_single_jam(jnr,type);% single freq
        Y = f_cal_fft(sj);
        b3(2,i,k) = f_par_b3(Y);
        nj = f_narrowbandnoise_jam(jnr,type); % narrow band noise jam
        Y = f_cal_fft(nj);
        b3(3,i,k) = f_par_b3(Y);        
        bj = f_broadbandnoise_jam(jnr,type);% broad band noise jam
        Y = f_cal_fft(bj);
        b3(4,i,k) = f_par_b3(Y);        
        cj = f_combnoise_jam(jnr,type);% comb-like noise jam
        Y = f_cal_fft(cj);
        b3(5,i,k) = f_par_b3(Y);     
        sj = f_scanning_jam(jnr,type);% comb-like noise jam
        Y = f_cal_fft(sj);
        b3(6,i,k) = f_par_b3(Y);   
    end
end
for  i = 1:6
    for j = 1:length(JNR)
        b3a(i,j) = mean(b3(i,j,:));
    end
end


row = size(b3a,1);
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:row
    plot(JNR,b3a(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('多音干扰','单音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('b3')
ylim([0,10])

% Chracater R in freq domain reflect the changing trend of waveform
function b3 = f_par_b3(Y)
% en = Y.^2;
en = Y;
uf = mean(en);
simgaf = sqrt(var(en));
b3 = mean((Y-uf).^3)/simgaf.^3;
end

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

四、 决策树法分类

基于决策树理论的分类是利用信号的特征参数与预设的门限$T$进行比较，大于门限的归为类别集合$A$，小于门限的归于集合$B$，不断进行分支，最终将信号归为某一类别。影响决策树性能的是判决的逻辑，一方面要选取合适的特征参数，一方面要合理确定阈值。
决策策略的选择，需要以下规则：1）不同JNR下，某一干扰信号对应的某一特征参数需要与另一类干扰信号的特征参数存在明显的差异，并且，为了保证分类性能，最好是满足$\min [{\mathbf{P}}_A]>\max [{\mathbf{P}}_B]$ , ${\mathbf{P}}_{\mathbf{A}}$、${\mathbf{P}}_{\mathbf{B}}$分别表示干扰信号$A$、$B$的某一参数在不同JNR下的向量。2）决策过程应该按照参数对于不同信号的差异程度进行，对于差异最明显的参数，应优先判决。
根据第三部分特征参数提取中的讨论，在分类中必须要选择能够对不同信号进行有效区分的特征参数，这里使用频域中的4个参数：载波因子系数C、平均频谱平坦系数Fse、Rf参数、频域矩偏度系数b3。
根据特征参数在不同JNR上的趋势变化，可以选择不同的决策策略。这里提供2种策略。其中策略1的决策树判决流程图如Figure 18所示。
门限选定：T©=5.04，T(b31)=7.03，T(Fse)=0.013，T(Rf)=2，T(b32)=4.2。
为方便实验，将干扰信号类型进行编号如下表：

采用策略1进行判决，得到混淆矩阵如Figure 19所示以及在不同JNR下的正确识别率如Figure 20所示。可以看出，在JNR在0~15dB范围内，对于除多音干扰外的全部干扰信号，正确识别率均大于95%。其中，单音干扰、宽带干扰、扫频干扰的正确识别率在不同的JNR下均为100%，梳状谱干扰的平均正确识别率达到99.9%。在JNR>2dB以后， 多音干扰正确识别率大于95%。根据混淆矩阵，造成多音干扰在低JNR时识别率低的问题主要是在低JNR时，由于b3参数比较接近，窄带干扰信号被识别为多音干扰信号。


策略2的决策树判决流程图如Figure 21所示。
门限选定：T©=5.04，T(Fse1)=0.013，T(Rf)=2，T(b3)=4.2，T(Fse2)=0.0215。

采用策略2进行判决，得到混淆矩阵如Figure 22所示以及在不同JNR下的正确识别率如Figure 22所示。可以看出，在JNR在0~15dB范围内，所有干扰信号的正确识别率均大于95%，性能较策略1有了明显提升，多音干扰与窄带干扰已经能够有效地被区分开。说明决策树分类的性能与决策策略有直接的关系，而决策策略又与特征参数的选取以及门限的确定有关。总之，决策树是一种需要人为决策的方法，在实现过程中容易由于人的主观因素，而无法达到最优分类。但是决策树的优点是，执行过程简单，复杂度低，适合特征参数差异比较大的情况。

% demo of judge on jam 
% analysis on the recognization rate
clear;close all;
JNR = 0:1:15;
mm = 100;% Mont-carlo times 求平均
typein = 1;
type = 2;%  1 静态生成 2 动态生成 
times = 1;
for j = 1:6
    typein = j;
    for i = 1:length(JNR)
        jnr = JNR(i);
        for k=1:mm                        
            tar_output(times) = j;
            [flag(j,i,k),output(times) ]= f_rec(typein,jnr,type);  
%             [flag(j,i,k),output(times) ]= f_rec2(typein,jnr,type);  
            times = times + 1;
        end
        flag_s(j,i) = sum(flag(j,i,:)/mm);        
    end    
end

tar_output_f = full(ind2vec(tar_output));
output_f = full(ind2vec(output));
plotconfusion(tar_output_f,output_f);

row = size(flag_s,1);
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:row
    plot(JNR,flag_s(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('单音干扰','多音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('正确识别率')
% ylim([0.88,1.02])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

function [flag,type] = f_rec(typein,jnr,type_par)
% type_par = 1 静态参数生成 type_par=2 动态生成
% typein = randi([1,6]);
% 输出 type 1 单音 2 多音 3 窄带 4 宽带 5 梳状 6 扫频
% 输出 flag = 1 正确 flag=0 错误识别
switch typein
    case 1
        x = f_single_jam(jnr,type_par); % test signal
    case 2
        x = f_multi_jam(jnr,type_par); % test signal
    case 3
        x = f_narrowbandnoise_jam(jnr,type_par); % test signal
    case 4
        x = f_broadbandnoise_jam(jnr,type_par); % test signal
    case 5
        x = f_combnoise_jam(jnr,type_par); % test signal
    case 6
        x = f_scanning_jam(jnr,type_par); % test signal
end
Y = f_cal_fft(x);% spectrum of test signal
% a4 = f_par_a4(x);
C = f_par_C(Y);
% A1= f_par_A1(Y);
Fse = f_par_Fse(Y);
Rf = f_par_R_f(Y);
b3 = f_par_b3(Y);

    if C>5.4
        type = 1;
    elseif b3>7.6
        type = 2;
    elseif Fse<0.013
    	type = 6;
    elseif Rf<2
        type = 4;
    elseif b3<4.2
        type = 5;
    else
        type = 3;
    end

% type_set = {'单音','多音','窄带','宽带','梳状','扫频'};

% input = type_set{typein}
% output = type_set{type}
if type-typein==0
    flag = 1; % 如果识别正确 返回1
else 
    flag = 0;
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

% 第2种策略
function [flag,type] = f_rec2(typein,jnr,type_par)
%
% type_par = 1 静态参数生成 type_par=2 动态生成
% typein = randi([1,6]);
% 输出 type 1 单音 2 多音 3 窄带 4 宽带 5 梳状 6 扫频
% 输出 flag = 1 正确 flag=0 错误识别
switch typein
    case 1
        x = f_single_jam(jnr,type_par); % test signal
    case 2
        x = f_multi_jam(jnr,type_par); % test signal
    case 3
        x = f_narrowbandnoise_jam(jnr,type_par); % test signal
    case 4
        x = f_broadbandnoise_jam(jnr,type_par); % test signal
    case 5
        x = f_combnoise_jam(jnr,type_par); % test signal
    case 6
        x = f_scanning_jam(jnr,type_par); % test signal
end
Y = f_cal_fft(x);% spectrum of test signal
% a4 = f_par_a4(x);
C = f_par_C(Y);
% A1= f_par_A1(Y);
Fse = f_par_Fse(Y);
Rf = f_par_R_f(Y);
b3 = f_par_b3(Y);

    if C>5.4
        type = 1;
    elseif Fse<0.013
        type = 6;
    elseif Rf<2
        type = 4;
    elseif b3<4.2
    	type = 5;    
    elseif Fse>0.0215
        type = 2;
    else
        type = 3;
    end

% type_set = {'单音','多音','窄带','宽带','梳状','扫频'};

% input = type_set{typein}
% output = type_set{type}
if type-typein==0
    flag = 1; % 如果识别正确 返回1
else 
    flag = 0;
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

五、 支持向量机（SVM）分类

支持向量机是一种典型的基于学习的分类方法。在有限个数据样本下，进行数据处理，得到已知训练数据下的最优解，然而将问题转化为一个二次寻优问题，得到全局最优解。
这里支持向量机的输入特征与决策树中一致，使用频域中的4个参数：载波因子系数C、平均频谱平坦系数Fse、Rf参数、频域矩偏度系数b3。核函数采用径向基核函数。

在每一干噪比下下，每一种干扰信号产生20个样本，作为训练数据。JNR范围为0~15dB，步进为1dB，对于每一干扰信号一共有320个样本。测试数据每一种干扰信号产生100个样本。使用Matlab自带的svmtrain函数进行训练，使用svmclassify函数进行分类。得到混淆矩阵如X所示以及在不同JNR下的正确识别率如Figure 25所示。可以看出，在JNR在0~15dB范围内，所有干扰信号的正确识别率均大于95%，性能较决策树法有了明显提升。但是相比决策树方法，SVM在分类中运算比较复杂，消耗时间比较长。

% type 1 单音 2 多音 3 窄带 4 宽带 5 梳状 6 扫频
%% 产生训练数据
num = 20;
JNR = 0:1:15;
len = 6*length(JNR)*num;
C_list = ones(len,1);Fse_list = ones(len,1);Rf_list = ones(len,1);
b3_list = ones(len,1);type_list = ones(len,1);
times = 1;
type = 2;%  1 静态生成 2 动态生成 
for typein =1:6
    for i = 1:length(JNR)
        jnr = JNR(i);
        for j = 1:num
            x = f_produce_jam(jnr,typein,type);% 产生干扰信号
            Y = f_cal_fft(x);% spectrum of test signal
            C_list(times) = f_par_C(Y);% 提取参数
            Fse_list(times) = f_par_Fse(Y);
            Rf_list(times) = f_par_R_f(Y);
            b3_list(times) = f_par_b3(Y);
            type_list(times) = typein;  % 数据标签，第x类干扰
            times = times+1;
        end
    end  
end
% 四组特征参数组成训练数据
tr_data = [C_list,Fse_list,Rf_list,b3_list,type_list];

%% 训练
tar_output = tr_data(:,5);%目标输出 数据标签
data = tr_data(:,1:4);% 特征参数
for i = 1:6
    Training2(i) = {data(find(tar_output==i),:)}; %将数据存入
end
num = nchoosek(1:6,2);% 6种类别两两组合
SVM = cell(size(num,1),1);%元胞形式的训练集
for k = 1:size(num,1)
    t1 = Training2{num(k,1)};
    t2 = Training2{num(k,2)};
    SVM{k} = svmtrain([t1;t2],[ones(size(t1,1),1);zeros(size(t2,1),1)],...
        'Kernel_function','rbf');%训练函数
end
disp('SVM训练完成')



%% 进行测试
JNR = 0:1:15;output_list=[];tar_output_list=[];
for i = 1:length(JNR)
    jnr = JNR(i);
    [output_temp,tar_output_temp,Pi(:,i)] = f_SVMtest(jnr,SVM,type);  % 6*16矩阵
    output_list = [output_list;output_temp];
    tar_output_list = [tar_output_list;tar_output_temp];
end

tar_output_f = full(ind2vec(tar_output_list'));
output_f = full(ind2vec(output_list'));
plotconfusion(tar_output_f,output_f);

row = size(Pi,1);
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:row
    plot(JNR,Pi(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('单音干扰','多音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('正确识别率')
xlim([0,15]);ylim([0.88 1.02])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

% 生成某一信噪比下的6种信号的测试数据 并进行测试
% 返回一个6*1的向量
function [output,type_list,Pi] = f_SVMtest(jnr,SVM,type)
num = 100;
len = 6*num;
C_list = ones(len,1);Fse_list = ones(len,1);Rf_list = ones(len,1);
b3_list = ones(len,1);type_list = ones(len,1);
times = 1;
for typein =1:6
        for j = 1:num
            x = f_produce_jam(jnr,typein,type);% 产生干扰信号
            Y = f_cal_fft(x);% spectrum of test signal
            C_list(times) = f_par_C(Y);% 提取参数
            Fse_list(times) = f_par_Fse(Y);
            Rf_list(times) = f_par_R_f(Y);
            b3_list(times) = f_par_b3(Y);
            type_list(times) = typein;  % 数据标签，第x类干扰
            times = times+1;
         end  
end
% 四组特征参数组成测试数据
te_data = [C_list,Fse_list,Rf_list,b3_list,type_list];
tar_output_test = te_data(:,5);% 目标输出结果
test = te_data(:,1:4);% 特征参数
num = nchoosek(1:6,2);% 6种类别两两组合
for kk = 1:size(test,1)
    for k = 1:length(SVM)
        result(k) = svmclassify(SVM{k},test(kk,:)); %分类
        temp(k) = num(k,1).*result(k)+num(k,2).*~result(k);
    end
    results(kk) = mode(temp,2);
end
output = results';

Pi = zeros(6,1);
for i = 1:6
    Posi = find(tar_output_test==i);
    deltai = tar_output_test(Posi)-output(Posi);
    Pi(i) =length(find(deltai==0))/size(Posi,1);
    % 统计不同干扰信号...的正确识别率
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

六、 神经网络（NN）分类

神经网络方法也是一种基于学习的算法，在模式识别
领域，通过训练得到的网络可以有效解决非线性分类问题。这里采用LM神经网络。识别结果如Figure 27及Figure 28所示。可以看出，将特征参数输入神经网络后进行分类，可以达到很好的识别率。在不同的JNR下，各类干扰信号均能实现100%的识别率。

% type 1 单音 2 多音 3 窄带 4 宽带 5 梳状 6 扫频
% 产生训练数据和测试数据
%% 产生训练数据
num = 20;
JNR = 0:1:15;
len = 6*length(JNR)*num;
C_list = ones(len,1);Fse_list = ones(len,1);Rf_list = ones(len,1);
b3_list = ones(len,1);type_list = ones(len,1);
times = 1;
type = 2;%  1 静态生成 2 动态生成 
for typein =1:6
    for i = 1:length(JNR)
        jnr = JNR(i);
        for j = 1:num
            x = f_produce_jam(jnr,typein,type);% 产生干扰信号
            Y = f_cal_fft(x);% spectrum of test signal
            C_list(times) = f_par_C(Y);% 提取参数
            Fse_list(times) = f_par_Fse(Y);
            Rf_list(times) = f_par_R_f(Y);
            b3_list(times) = f_par_b3(Y);
            type_list(times) = typein;  % 数据标签，第x类干扰
            times = times+1;
        end
    end  
end
% 四组特征参数组成训练数据
% 格式： 1:4列为特征参数 第5列为数据标签 表示第n类干扰信号
tr_data = [C_list,Fse_list,Rf_list,b3_list,type_list];
save('training_data.mat','tr_data');
disp('训练数据产生完成')

%% 生成测试数据
JNR = 0:1:15;
num = 20;% 每一个JNR产生num个测试数据
len = 6*length(JNR)*num;%6种干扰信号
C_list = ones(len,1);Fse_list = ones(len,1);Rf_list = ones(len,1);
b3_list = ones(len,1);type_list = ones(len,1);JNR_list = ones(len,1);
times = 1;    
for typein =1:6
    for i = 1:length(JNR)
        jnr = JNR(i);
        for j = 1:num  % 每一个JNR下有20个数据
            x = f_produce_jam(jnr,typein,type);% 产生干扰信号
            Y = f_cal_fft(x);% spectrum of test signal
            C_list(times) = f_par_C(Y);% 提取参数
            Fse_list(times) = f_par_Fse(Y);
            Rf_list(times) = f_par_R_f(Y);
            b3_list(times) = f_par_b3(Y);
            type_list(times) = typein;  % 数据标签，第x类干扰
            JNR_list(times) = jnr;  % 数据标签，JNR                
            times = times+1;
        end  
    end
end
te_data = [C_list,Fse_list,Rf_list,b3_list,JNR_list,type_list];

save('test_data.mat','te_data');
disp('测试数据产生完成')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

% type 1 单音 2 多音 3 窄带 4 宽带 5 梳状 6 扫频
% 神经网络训练和测试
%% 训练
clear all;close all;
load('training_data.mat'); % 导入训练数据
disp('训练数据已导入')
input = tr_data(:,1:4);% 属性数据
tar_output = tr_data(:,5);% 标签数据

% 数据变换
input = input';
tar_output = tar_output';
tar_output = full(ind2vec(tar_output));

% 设置LM网络参数
net = patternnet(10,'trainlm');
net.trainParam.epochs = 1000;
net.trainParam.show = 25;
net.trainParam.showCommandLine = 0;
net.trainParam.showWindow = 0;
net.trainParam.goal = 0;
net.trainParam.time = inf;
net.trainParam.min_grad = 1e-6;
net.trainParam.max_fail = 5;
net.performFcn = 'mse';

net = train(net,input,tar_output);
output = sim(net,input);
figure;plotconfusion(tar_output,output);% 混淆图

save('net.mat','net');
disp('训练完成')

output = vec2ind(output);
output=output';
ori_label = tr_data(:,5);% 标签数据

flag = zeros(length(output),1);
for i = 1:length(output)
    if output(i)-ori_label(i)==0
        flag(i)=1;
    else 
        flag(i)=0;
    end
end

for i = 1:6
    for j = 1:16
        left = (i-1)*320+20*j-19;
        right = (i-1)*320+20*j;
        rate(i,j) = sum(flag(left:right))/20;
    end
end

JNR = 0:1:15;
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:6
    plot(JNR,rate(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('单音干扰','多音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('训练数据正确识别率')
xlim([0,15]);ylim([0.88 1.02])

%% 测试
load('test_data.mat'); % 导入训练数据
load('net.mat'); % 导入训练数据
disp('测试数据已导入')
input = te_data(:,1:4);% 属性数据
tar_output = te_data(:,6);% 标签数据

output = sim(net,input');
tar_output = full(ind2vec(tar_output'));
figure;plotconfusion(tar_output,output);% 混淆图

output = vec2ind(output);
output=output';
ori_label = te_data(:,6);% 标签数据

flag = zeros(length(output),1);
for i = 1:length(output)
    if output(i)-ori_label(i)==0
        flag(i)=1;
    else 
        flag(i)=0;
    end
end

for i = 1:6
    for j = 1:16
        left = (i-1)*320+20*j-19;
        right = (i-1)*320+20*j;
        rate(i,j) = sum(flag(left:right))/20;
    end
end

JNR = 0:1:15;
figure;lintype = {'-b^','-bv','-b<','-b>','-bs','-bo'};
colortype = {[1,0.1,0.1],[1,0.7,0],[1,1,0],[0,0.5,0],...
    [0,0.8,0.8],[1,0.08,1]};
for i = 1:6
    plot(JNR,rate(i,:),char(lintype(i)),'MarkerFaceColor',colortype{i});
    hold on;
end
hold off;grid on;
legend('单音干扰','多音干扰','窄带干扰','宽带干扰','梳状谱干扰','扫频干扰')
xlabel('JNR(dB)')
title('测试数据正确识别率')
xlim([0,15]);ylim([0.88 1.02])

disp('测试完成')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

补充代码

% Chracater the spectrum of jam x 
function Y = f_cal_fft(x)
         % Sampling frequency
x_mid = (max(x)+min(x))/2;
x_nor = 2*(x-x_mid)/(max(x)-min(x));% normalize the data into [-1,1]
x_dec = (x_nor-mean(x_nor))/sqrt(var(x_nor));% decentralize the data
x = x_dec;
L = length(x);% signal length
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(x,NFFT)/L;
Y = 2*abs(Y(1:NFFT/2+1));%single-wise of spectrum


1
2
3
4
5
6
7
8
9
10
11
12
13