深 圳 大 学 实 验 报 告
实验课程名称： 人工神经网络技术
实验项目名称： BP神经网络对蝴蝶花分类
学院： 专业： 软件工程
报告人： 学号： 班级：
同组人： 无
指导教师：
实验时间：
实验报告提交时间：
教务处制
一、 实验目的
初步熟悉MATLAB 工作环境，熟悉命令窗口，学会使用帮助窗口查找帮助信息。
二、 实验内容
1、网络设计，包括输入层、隐含层、输出层节点个数的设计。
2、算法步骤
3、编程，注意原始数据的通用化，数据输入的随机性。
4、网络训练，注意训练数据与验证数据分开。
5、网络验证
6、结果分析，修改隐含层节点个数，修改学习率，分别对结果的影响。
三、实验步骤
直接在Matlab软件中的Editor中运行以下代码：（完善的代码用红色字体表示）
% li.m
% A BP NN with momentum to solve Fisher's Iris Flower problem
% by lixiujuan, Nov 13, 2011
%
% the NN architecture
% it is a three layers neural network 4-3-3.
%
% parameter description
% h=4 the node numer of input layer
% i=3 the node numer of hidden layer
% j=3 the node numer of output layer
% V[h,i] the weights between input and hidden layers
% W[i,j] the weights between hidden and output layers
% Pi[i] the thresholds of hidden layer nodes
% Tau[j] the thresholds of output layer nodes
% a[h] the input values
% b[i] the hidden layer node activations
% c[j] the output layer node activations
% ck[j] the desired output of output layer nodes
% d[j] the eror in output layer nodes
% e[i] the eror in hidden layer nodes
% DeltaW[i,j] the amount of change for the weights W[i,j]
% DeltaV[h,i] the amount of change for the weights V[h,i]
% DeltaPi[i] the amount of change for the thresholds Pi[i]
% DeltaTau[j] the amount of change for the thresholds Tau[j]
% Alpha=0.1 the leaning rate
% Beta=0.1 the leaning rate
% Gamma=0.8 the constant determines effect of past weight changes
% Tor=0.001 the torrelance that determines when to stop training
% Maxepoch=1000 the max iterate number
%
% other parameters
% Ntrain=115 the number of trainning sets
% Ntest=35 the number of test sets
% Otrain[115] the output of training sets
% Otest[35] the output of test sets
% Odesired[150] the desired output of training and test sets
%
% function description
% f(x)=logsig(x)
% f(x)=1/(1+exp(-x))
%
% data file
% input file: data.dat
%
close all; clc; clf; clear all;
% parameters for the NN structure
h=4;
i=3;
j=3;
Alpha=0.1;
Beta=0.1;
Gamma=0.85;
Tor=0.0005;
Maxepoch=2000;
Accuracy=0;
Ntrain=115;
Ntest=35;
%assign random values in the range [-1, +1]
V=2*(rand(h,i)-0.5);
W=2*(rand(i,j)-0.5);
Pi=2*(rand(1,i)-0.5);
Tau=2*(rand(1,j)-0.5);
DeltaWOld(i,j)=0; %set the delat of Wij to 0
DeltaVOld(h,i)=0; %set the delat of Vij to 0
DeltaPiOld(i)=0; %set the delat of Pi to 0
DeltaTauOld(j)=0; %set the delat of Tau to 0
% the learning process
Epoch=1;
Error=10;
% load the training set data and test set data
load data.dat
Odesired=data(:,2); % get the desired of output of 150 data sets
% normalize the input data to rang [-1 +1]
datanew=data(:,3:6);
maxv=max(max(datanew));
minv=min(min(datanew));
datanorm=2*((datanew-minv)/(maxv-minv)-0.5);
while Error＞Tor
Err(Epoch)=0;
for k=1:Ntrain % k = the index of tranning set
a=datanorm(k,:); % get the input
% set the desired output ck[j]
if data(k,2)==0
ck=[1 0 0];
elseif data(k,2)==1
ck=[0 1 0];
else
ck=[0 0 1];
end;
% calculate the hidden nodes activation
for ki=1:i
b(ki)=logsig(a*V(:,ki)+Pi(ki));
end;
% calculate the output nodes activation
for kj=1:j
c(kj)=logsig(b*W(:,kj)+Tau(kj));
end;
% calculate error in output Layer FC
d=c.*(1-c).*(ck-c);
% calculate error in hidden layer FB
e=b.*(1-b).*(d*W');
% adjust weights Wij between FB and FC
for ki=1:i
for kj=1:j
DeltaW(ki,kj)=Alpha*b(ki)*d(kj)+Gamma*DeltaWOld(ki,kj);
end
end;
W=W+DeltaW;
DeltaWOld=DeltaW;
% adjust weights Vij between FA and FB
for kh=1:h
for ki=1:i
DeltaV(kh,ki)=Beta*a(kh)*e(ki);
end
end;
V=V+DeltaV;
DeltaVold=DeltaV;
% adjust thresholds Pi and Tau
DeltaPi=Beta*e+Gamma*DeltaPiOld;
Pi=Pi+DeltaPi;
DeltaPiold=DeltaPi;
DeltaTau=Alpha*d+Gamma*DeltaTauOld;
Tau=Tau+DeltaTau;
DeltaTauold=DeltaTau;
% the error is the max of d(1),d(2),d(3)
Err(Epoch)=Err(Epoch)+0.5*(d(1)*d(1)+d(2)*d(2)+d(3)*d(3));
end %for k=1:Ntrain
Err(Epoch)=Err(Epoch)/Ntrain;
Error=Err(Epoch);
% the training stops when iterate is too much
if Epoch ＞ Maxepoch
break;
end
Epoch = Epoch +1; % update the iterate number
end
% test data
for k=1:Ntest % k = the index of test set
a=datanorm(Ntrain+k,:); % get the input of test sets
% calculate the hidden nodes activation
for ki=1:i
b(ki)=logsig(a*V(:,ki)+Pi(ki));
end;
% calculate the output of test sets
for kj=1:j
c(kj)=logsig(b*W(:,kj)+Tau(kj));
end;
% transfer the output to one field format
if (c(1)＞ 0.9)
Otest(k)=0;
elseif (c(2)＞ 0.9)
Otest(k)=1;
elseif (c(3)＞ 0.9)
Otest(k)=2;
else
Otest(k)=3;
end;
% calculate the accuracy of test sets
if Otest(k)==Odesired(Ntrain+k)
Accuracy=Accuracy+1;
end;
end; % k=1:Ntest
% plot the error
plot(Err);
% plot the NN output and desired output during test
N=1:Ntest;
figure; plot(N,Otest,'b-',N,Odesired(116:150),'r-');
% display the accuracy
Accuracy = 100*Accuracy/Ntest;
t=['TESTING RESULT, the accuracy of test sets is: ' num2str(Accuracy) '%' ];
disp(t);
当中间隐藏层结点数为i=3时，得出下图：
当Alpha=0.08;Beta=0.08时;如下图：
当Alpha=0.05;Beta=0.05时;如下图：
当Alpha=0.2;Beta=0.2时;如下图：
当Alpha=0.3;Beta=0.3时;如下图：
当Alpha=0.2;Beta=0.3时;如下图：
当Alpha=0.3;Beta=0.2时;如下图：
当中间隐藏层结点数为i=4时，得出下图：
当中间隐藏层结点数为i=2时，得出下图：
当Alpha=0.05;Beta=0.05时;如下图：
当中间隐藏层结点数为i=1时，得出下图：
四、总结分析
由于对MATLAB的不熟悉，一开始还是不知道怎么完成这个实验，但是在同学的帮助下，渐渐懂得怎么导入数据和代码实现，进而了解怎么样使用BP神经网络来进行蝴蝶花的分类。通过修改中间隐藏层的结点数，进行对比，发现i=3以及i=4的时候所得到的图是一样的，说明在i大于等于3的时候，实验结果误差较小以至于可以忽略不计，又进行了一次i=2的实验，结果发现结果也是几乎一样的，再进行进一步的实验，将i改为1，发现期望输出与实际输出存在明显区别，到后期才重合。所以，我觉得在这次的实验中，中间隐藏层结点数需要大于等于2比较合适。还有就是学习率的问题，刚开始实验，就只是改了隐藏节点，而忽略了学习率的问题。所以在原来的基础上，有改动了一下学习率，发现Figure1还是有很大地变动的，但是Figure2反而没什么不同。我觉得学习率的范围应该在0.1以内的。可是还是不知道学习率具体是起什么作用的。这次实验还是大概的完成了，但是对MATLAB还是有一些不懂的地方，仍需在接下来的实验中继续学习。
第二篇：实验报告
电子情报侦察分析综合实验
一、实验目的
领会专业的特色，掌握电子情报侦察分析的基本流程，理解雷达原理与雷达截获分析课程的联系。
二、实验要求
1、工具采用MATLAB
2、对实验结果进行理论分析，并给出自己的见解
三、实验计划
1、拟将18名学员分成两组，进行电子情报模拟对抗
2、每名学员编程产生雷达信号，汇集给教员，由教员随机配对，对产生的信号进行分析。
四、实验步骤和内容
4.1产生不同体制的雷达信号
1、实验原理
加强对简单脉冲信号、频率分集信号、频率捷变信号、重频参差信号、PRI抖动和PRI跳变信号、PRI滑变信号、PRI排定信号、脉组PRI变化信号、双脉冲信号、脉冲压缩信号等信号的直观认识；使学员从多角度描述雷达信号。
常见的雷达信号形式内容：
（1）简单脉冲信号
指载频、脉冲重复频率和脉冲宽度等三个参数均固定不变的信号
（2）频率分集信号
指一部雷达同时或在相隔很短的时间内发射多个载频以完成同一任务。频率分集信号中，不同载频的脉冲可以同时发射，也可以顺序发射。频率分集信号的频谱是各个不同频率脉冲频谱的合成。
（3）脉冲重复间隔（PRI）变化信号
PRI变化信号脉冲指脉冲的PRI在相邻脉冲之间或脉组之间变化的信号。
①重频参差信号。重频参差信号指PRI在相邻脉冲之间或脉组之间变化的信号
②PRI抖动信号和PRI跳变信号
PRI抖动指相邻脉冲的PRI在一定范围内抖动，即相邻的脉冲时间间隔不相等。一般PRI抖动的范围小于中心值的5%。
PRI跳变信号是指PRI在其平均值附近的有限几个点上按照某一分布规律跳动，其跳动范围可高达平均值得30%。
PRI抖动信号与PRI跳变信号的差别是，前者可以在无限个点上变化，后者只能在有限个点上变化；前者变化范围小，后者变化范围大。
③PRI滑变信号
这种信号的特点是PRI值周期性单调增加或下降，在两个极值之间周期性变化，PRI滑变信号可以消除盲距，还可以提供恒定高度范围内的仰角扫描，达到最佳性能。在这种应用中，最小PRI的值与最大的值之比大致等于仰角扫描系统的最小作用距离与最大作用距离之比。通常这个值为3~6。
④PRI排定信号
这种信号的特点是PRI值按排定的多个数值周期性的变化，其PRI数值一把为十几至二十几点。PRI排定信号一把用于多功能雷达。
⑤脉组PRI变化
脉组PRI变化信号是指信号在PRI为T1值上发射几个脉冲，然后又在PRI为T2值上发射几个脉冲，并且具有周期性。每组的脉冲数一般不少于四个。
（4） 双脉冲信号
这是指在每个PRI内发射两个脉冲，这两个脉冲的宽度可以相同，也可以不同，两脉冲间的间隔通常为脉冲宽度的几倍。
（5） 脉冲压缩信号
是指单脉冲内的频率呈线性变化。
2、实验要求
（1）、绘制出各种形式的雷达信号；
（2）、频率变化的信号要频率分析各雷达信号；
（3）、总结该雷达信号特征，说明选用该波形的目的以及这样波形的优缺点；
（4）、给出源代码。
4.2雷达信号截获与定位实验
1、实验原理：
不同雷达发射不同体制雷达信号，经电磁空间传播被雷达信号截获接收机接收，本实验，根据侦察方程模拟该过程。
仿真3~5部雷达信号，体制分别为简单体制、重频参差、脉压信号（线性调频）、频率捷变等，脉冲个数至少有100个；雷达信号统一为1；雷达天线增益函数为
，设每部雷达波束20度，扫描速度10转/分；不考虑损耗与大气传播影响；
雷达信号截获接收机有1~3个，且同时截获雷达信号；
截获接收机天线采用全向天线；
2、实验要求：
（1）、仿真出雷达信号截获接收机截获的雷达信号；
（2）、给出雷达信号参数（载频、脉宽、到达时刻、脉幅）测量方法，并给出测量值；
（3）、根据截获的雷达信号，给出相应的定位方法，并给出位置信息。
4.3雷达信号分选实验
1、实验原理
对截获接收机的PDW输出流进行去交错处理，即雷达信号脉冲分选。按照分选采用的参数的不同分为多参数分选与PRI分选等。
多参数分选即按照到达方向、载频、脉宽等参数进行分选，也可以看作是多参数聚类分析。常用的聚类算法有K均值、最大最小方法以及神经网络方法等。
PRI分选即在多参数分选后，按照雷达脉冲到达时刻信息，获取脉冲间隔信息，然后采用统计、变换域等方法进行PRI分选。
2、实验内容
（1）选取任一种聚类方法，完成多参数雷达信号多参数分选；
（2）采用CDIF等PRI分选方法，完成PRI分选，重点观察重频参差信号的分选效果。
3、实验要求
（1）以实验二的结果进行该实验。
五、实验结果及分析
5.1产生不同体制的雷达信号
1.简单脉冲信号
产生的简单脉冲如图5.1所示。
图5.1
参数设置如下：
2.线性调频脉冲信号
产生的脉冲如图5.2所示。
图5.2
参数设置如下：
3.重频参差脉冲信号
产生的脉冲如图5.3所示。
图5.3
参数设置如下：
4.重频参差脉冲信号
产生的脉冲如图5.4所示。
图5.4
参数设置如下：
1.仿真出雷达信号截获接收机截获的雷达信号；
（1）仿真两个雷达接收机，其参数如下：
一号接收机：（50km ,62km） 天线有效面积：1
二号接收机：（1.5km ,6km） 天线有效面积：1.2
在信号接收过程中，由于四部雷达的距离不同，会产生不同的时延，且在传输过程中会夹杂着噪声，因此，在仿真接收信号中加入了时延和高斯噪声（信噪比为20，信噪比不能太大，否则有可能噪声将信号淹没，无法进行检测）。
加时延的代码如下：
c=3e8;
shiyan1=rand(1,4)/1e3;
for i=1:4
shiyan2(i)=shiyan1(i)+(R2(i)-R1(i))/c;
end
加噪声的代码如下：
SNR=20;
snr=1/(10^(SNR/10));
Y11=[zeros(1,fs*shiyan1(1)) Y1];len11=length(Y11);
noise11=sqrt(snr/2)*(randn(1,len11));
Y_11=Y11+noise11;
在接收信号时，由于四种信号的点数不同，需要以最长的信号为准，在其他信号后补零，以保证正确的叠加。
（2）两个接收机收到的雷达信号分别如图5.5和图5.6所示。
图5.5
图5.6
2． 给出雷达信号参数（载频、脉宽、到达时刻、脉幅）测量方法，并给出测量值。
（1）从接收到的信号中提取脉冲
思路步骤：先设定一个门限值，然后对信号进行判别，如果大于门限，则判为起始位置，若小于门限，则判为结束位置。其中需要注意的是：
A．门限的设置必须合理；
B．在判断起始位置和结束位置时，需要考虑脉内小于或大于门限的点，因此，通过将后边200个点求和，设定一个门限，判断是否为起始和结束位置。
为了判断门限更简便，首先将信号二值化处理，然后比较容易的设定门限。
信号二值化代码如下：
as=zeros(1,length(s));
for i=1:length(s)
if s(i)＞=T %%如果信号幅度大于等于门限，则判为1
as(i)=1;
else
as(i)=0;%否则判为0
end
end
求脉冲的起始位置和结束位置的代码如下：
for i=2:length(s)-300
summ=sum(as(i:i+300));
if as(i)==1&&count==0
if summ＞=220
b(m,1)=i; % b矩阵的第一列存放脉冲起始位置
m=m+1;
count=1;
end
end
if count==1&&summ＜40
b(n,2)=i; % b矩阵的第二列存放脉冲结束位置
n=n+1;
count=0;
end
end
实验结果如图5.7所示：测得脉冲个数为212个
图5.7
（2）求脉宽
思路步骤：上一步将脉冲的起始位置和结束位置存放在一个矩阵b中，脉宽的测量只需用结束位置减去起始位置，求出点数，然后除以采样频率fs，即可得到脉宽。
求脉宽的代码如下：
a=b(:,2)-b(:,1); %a矩阵存放脉宽
figure;
stem(a);
axis([0,200,0,4000]);
实验结果如图：
图5.8
通过观察可看出，测得的脉宽大量集中于6us，10us，20us，30us，还有部分零散分布在30—40us，与了雷达发射信号的脉宽近似符合，测量结果比较准确。
（3）求载频
思路步骤：求载频时，将每个脉冲做傅里叶变换，从中提取出最大幅度的即为载频位置，然后根据公式fc=k/N*fs得到载频。
需要注意的是，在傅里叶变换后，会有较大的直流分量的影响，因此需要将直流分量去除，然后再求最大值。
求载频的代码如下：
N=10000;
for i=1:length(b)
maichong=[s(b(i,1):b(i,2))];
maichong_FFT=abs(fft(maichong,N));
maichong_FFT(1:60)=0;
maichong_FFT(9940:10000)=0;
[Fcmax,k]=max(maichong_FFT);
Fc(i)=k/N*fs; %Fc矩阵存放载频
End
实验结果如图5.9所示：
图5.9
通过观察可看出，测得的载频大量集中于5MHz，10MHz，15MHz，20MHz，还有部分零散分布在20—25MHz，与了雷达发射信号的载频近似符合，测量结果比较准确。
5.3雷达信号分选实验
思路步骤：根据测得的脉宽和载频，利用K均值聚类法，对测得的脉冲进行聚类分选，其难点在于K均值聚类。
K均值聚类的思想是，先人为选定四个与真实值较为接近的聚类中心，然后设定四个聚类种类，求出其余各点与聚类中心的距离，将这些距离求均值，如果此均值与原聚类中心的距离小于一个门限值，则说明此聚类中心比较合适，聚类结束；若此均值与原聚类中心的距离大于门限，则将此均值当做新的聚类中心，从头开始，在进行聚类，直到聚类结束，最终得到四类。
实验中较大的难点在于门限值的设置，需要根据每次得到的数据进行分析，设定一个合适的门限。
实验结果：
通过上表可以看出，四种体制的脉冲信号的脉宽和载频的测量值与实际值在误差允许的范围内较为精确，基本能够分选出四种体制信号的参数；但脉冲个数的聚类结果与实际值相差较大。
脉冲个数误差较大的原因主要有：a. 脉冲的重叠，产生很大脉宽的脉冲，在筛选脉冲时将其舍弃;
b. 聚类时，由于误差的存在，部分脉冲与聚类中心距离较大，在门限值之外，因此将其舍弃。
六、实验心得
本次实验是一个比较综合的实验，是对以前学习的五门课程知识的综合总结。虽然最终没有完成全部的实验内容，但通过完成的实验内容，我感觉收获已经很多。
实验中产生四种不同体制的雷达信号，这在电子战实验中已经做过了，但需要产生100个脉冲，因此需要理解脉冲重复间隔和时间的关系，根据采样定理，还需要使采样率大于2倍的载频。在这部分中，我对各种体制雷达信号的原理有了一定掌握，之前对于这些信号调制原理根本不懂。此外，在测量脉宽和载频时，加深了对信号点数与脉宽、载频之间的转换关系。
通过实验，我收获最大的是使用matlab编程的能力有了明显提高，这个实验的程序应该是我大学四年编程最长最综合的一次，在编程过程中不断理解原理，思考编程思路，如何将原理转变成程序来实现，不断调试程序、修改程序，最终得到较为满意的实验结果，这为我下学期的毕业设计打下了较好的基础。
实验过程中，我得到了教员以及一些能力水平比较高的同志的帮助和指导，可以看出自己与某些同志之间的差距，这对我的帮助也是很大的，在感谢教员和同志们的同时，我也会从整个实验过程中总结经验，在今后的学习中弥补不足。
七、源代码附录
%% %%%%%%%%简单脉冲雷达信号%%%%%%%%
clc
clear all
close all
PW=6e-6; %脉宽
fs=10e7; %采样频率
fc=10e6; %载频
prt=30e-6;%脉冲重复间隔
pri=1/prt; %脉冲重复频率
t=0:1/fs:100*prt;
y=cos(2*pi*fc*t-0.5);
duty=100*PW/prt; %占空比
s=(square(2*pi*pri*t,duty)+1)/2; %产生方波
Y1=y.*s;
plot(t*1e6,Y1);
xlabel('t/us');ylabel('幅度');title('简单脉冲信号');
%% %%%%%%%%%%%线性调频雷达信号%%%%%%%%%%
pw=20e-6; %发射脉宽
B=10e6; %调频带宽
% fs=10e7; %采样频率
fc=20e6; %载频
prt=70e-6;
pri=1/prt; %脉冲重频
t=0:1/fs:prt;
f0=rectpuls(t-pw/2,pw);
Fc=f0.*exp(j*2*pi*(fc*t+B/pw*t.*t));
% figure
% plot(abs(fft(Fc)));
% figure;
% plot(t*1e6,real(Fc));title('线性调频信号');xlabel('t/us');
Y2=repmat(Fc,1,100);
figure;
plot(real(Y2));title('线性调频脉冲发射信号');
%% %%%%%%%%%PRI参差%%%%%%%%%
fc=5e6; %载频
tao=10e-6; %脉宽
prt1=100e-6;
prt2=200e-6;
prt3=300e-6; %脉冲间隔
% fs=3*fc; %采样频率
prt=prt1+prt2+prt3; %参差脉冲间隔
pri=1/prt; %参差脉冲重频
t=0:1/fs:33*prt;
zhankong=tao/prt*100; %占空比
x=sin(2*pi*fc*t);
y =(square(2*pi*pri*t,zhankong)+1)/2;
s1=(square(2*pi* pri *(t-prt1), zhankong)+1)/2;
s2=(square(2*pi* pri *(t-prt1-prt2), zhankong)+1)/2;
s3=(square(2*pi* pri *(t-prt1-prt2-prt3), zhankong)+1)/2;
Y3=x.*s1+x.*s2+x.*s3;
figure;
plot(t*1e6,abs(Y3));title('PRI参差信号');
xlabel('t/us');ylabel('幅度');
%% % %%%%%%%重频抖动信号%%%%%%
fc=15e6; %载频
pw=30e-6;%脉宽
prix=(rand(1,4)*0.5*50+50)*1e-6;
pri1=prix(1);pri2=prix(2);pri3=prix(3);pri4=prix(4);
pri=pri1+pri2+pri3+pri4;
prt=1/pri;
zhankong=pw/pri*100; %占空比
t=0:1/fs:25*pri;
x=sin(2*pi*fc*t);
y1=(square(2*pi*prt*t,zhankong)+1)/2;
y2=(square(2*pi*prt*(t-pri1),zhankong)+1)/2;
y3=(square(2*pi*prt*(t-pri1-pri2),zhankong)+1)/2;
y4=(square(2*pi*prt*(t-pri1-pri2-pri3),zhankong)+1)/2;
y5=(square(2*pi*prt*(t-pri1-pri2-pri3-pri4),zhankong)+1)/2;
Y4=x.*y2+x.*y3+x.*y4+x.*y5;
figure;
plot(t*1e6,abs(Y4));
title('PRI抖动雷达脉冲信号');
xlabel('t/us');ylabel('幅度');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 雷达信号截获
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%4
%%%设定雷达和接收机的坐标
% clc;close all;clear all;
leidax=[1e4 2.3e4 3.5e4 7e4];
leiday=[4e4 1.8e4 4.5e4 3e4];
jieshoux=[5e4 6.2e4];
jieshouy=[1.5e3 6e3];
% stem(leidax,leiday);hold on;
% stem(jieshoux,jieshouy,'r');
% axis([0,1e5,0,5e4]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 计算接收机与雷达的距离
%%%%%%%%与接收机1的距离%%%%%%
for i=1:4
R1(i)=sqrt((leidax(i)-jieshoux(1))^2+(leiday(i)-jieshouy(1))^2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%与接收机2的距离%%%%%%
for i=1:4
R2(i)=sqrt((leidax(i)-jieshoux(2))^2+(leiday(i)-jieshouy(2))^2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pt=[1e7 2e7 3e7 4e7];%各雷达的发射功率
Ar=[1 1.2];%两个接收机的有效接收面积
Gt=[30 31 32 33];
Gtt=10.^(Gt./10); %各雷达天线增益
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 接收机的接收功率
for i=1:4
Pr1(i)=(Pt(i).*Gtt(i)*Ar(1))./(4*pi*(R1(i)^2));%求接收机1的接收功率
end
for i=1:4
Pr2(i)=(Pt(i).*Gtt(i)*Ar(2))./(4*pi*(R1(i)^2));%求接收机2的接收功率
end
%% 加时沿
c=3e8;
shiyan1=rand(1,4)/1e3; %产生接收机1的随机时延
for i=1:4
shiyan2(i)=shiyan1(i)+(R2(i)-R1(i))/c; %接收机2的时延与时延1有关
end
%% 加噪声
SNR=20; %信噪比
snr=1/(10^(SNR/10));
% noise=-30;
%%%%接收机1接收的信号
Y11=[zeros(1,fs*shiyan1(1)) Y1];len11=length(Y11); %%加时延
Y12=[zeros(1,fs*shiyan1(2)) Y2];len12=length(Y12);
Y13=[zeros(1,fs*shiyan1(3)) Y3];len13=length(Y13);
Y14=[zeros(1,fs*shiyan1(4)) Y4];len14=length(Y14);
% noise11=awgn(Y11,noise);
% noise12=awgn(Y12,noise);
% noise13=awgn(Y13,noise);
% noise14=awgn(Y14,noise);
noise11=sqrt(snr/2)*(randn(1,len11));
noise12=sqrt(snr/2)*(randn(1,len12));
noise13=sqrt(snr/2)*(randn(1,len13));
noise14=sqrt(snr/2)*(randn(1,len14));
Y_11=Y11+noise11; %%加噪声
Y_12=Y12+noise12;
Y_13=Y13+noise13;
Y_14=Y14+noise14;
YY1=[Y_11 zeros(1,length(Y_13)-length(Y_11))];
YY2=[Y_12 zeros(1,length(Y_13)-length(Y_12))];
YY4=[Y_14 zeros(1,length(Y_13)-length(Y_14))];
YY3=Y_13;
receive1=YY1+YY2+YY3+YY4; %%接收机1接收的信号
figure;
plot(abs(receive1));
title('接收机1收到的信号');
%%%%接收机2接收的信号
Y21=[zeros(1,fs*shiyan2(1)) Y1];len21=length(Y21); %加时延
Y22=[zeros(1,fs*shiyan2(2)) Y2];len22=length(Y22);
Y23=[zeros(1,fs*shiyan2(3)) Y3];len23=length(Y23);
Y24=[zeros(1,fs*shiyan2(4)) Y4];len24=length(Y24);
% noise21=awgn(Y21,noise);
% noise22=awgn(Y22,noise);
% noise23=awgn(Y23,noise);
% noise24=awgn(Y24,noise);
noise21=sqrt(snr/2)*(randn(1,len21));
noise22=sqrt(snr/2)*(randn(1,len22));
noise23=sqrt(snr/2)*(randn(1,len23));
noise24=sqrt(snr/2)*(randn(1,len24));
Y_21=Y21+noise21;
Y_22=Y22+noise22; %加噪声
Y_23=Y23+noise23;
Y_24=Y24+noise24;
YYY1=[Y_21 zeros(1,length(Y_23)-length(Y_21))];
YYY2=[Y_22 zeros(1,length(Y_23)-length(Y_22))];
YYY4=[Y_24 zeros(1,length(Y_23)-length(Y_24))];
YYY3=Y_23;
receive2=YYY1+YYY2+YYY3+YYY4; %%接收机2接收的信号
figure;
plot(abs(receive2));
title('接收机2收到的信号');
fid = fopen('C:\Documents and Settings\Administrator\桌面\信号分析\receive1.dat','w');
fwrite(fid,receive1,'double');
fid = fopen('C:\Documents and Settings\Administrator\桌面\信号分析\receive2.dat','w');
fwrite(fid,receive2,'double');
%% 信号分析
clc;close all;clear all;
fid=fopen('C:\Documents and Settings\Administrator\桌面\信号分析\receive1.dat');
s=fread(fid,'double');
s=abs(s');
plot(s);
% noise=s(1:10000);
% noisemean=mean(noise);
% noisevar=var(noise);
T=0.3; %%%设定门限
fs=10e7;
% b=zeros(length(s),2);
%%脉宽
m=1;n=1;q=0;summ=0;count=0;%定义参数初值
%% 将信号二值化
as=zeros(1,length(s));
for i=1:length(s)
if s(i)＞=T %%如果信号幅度大于等于门限，则判为1
as(i)=1;
else
as(i)=0;%否则判为0
end
end
%% 求脉冲的起始位置和结束位置
for i=2:length(s)-300
summ=sum(as(i:i+300));
if as(i)==1&&count==0
if summ＞=220
b(m,1)=i; % b矩阵的第一列存放脉冲起始位置
m=m+1;
count=1;
end
end
if count==1&&summ＜40
b(n,2)=i; % b矩阵的第二列存放脉冲结束位置
n=n+1;
count=0;
end
end
%% 求脉宽
a=b(:,2)-b(:,1); %a矩阵存放脉宽
figure;
stem(a);
axis([0,200,0,4000]);
title('脉宽');
%% 载频
N=10000;
for i=1:length(b)
maichong=[s(b(i,1):b(i,2))];
maichong_FFT=abs(fft(maichong,N));
maichong_FFT(1:60)=0;
maichong_FFT(9940:10000)=0;
[Fcmax,k]=max(maichong_FFT);
Fc(i)=k/(2*N)*fs; %Fc矩阵存放载频
end
figure;
stem(Fc);
title('载频');
%% K均值聚类法对脉冲进行分类
% %选择聚类中心
maizu=[a';Fc/10000];
jlzx=[maizu(:,1),maizu(:,2),maizu(:,102),maizu(:,128)]; %聚类中心
jlzl(1,1)=1;
jlzl(2,1)=2;
jlzl(3,1)=102;
jlzl(4,1)=128; %聚类种类
%%
%设置判决条件
panjue=20;
panjue_js=25;
while panjue_js＞panjue
%对每个脉冲进行分类
jljl=zeros(1,4); %聚类距离
m1=2;m2=2;m3=2;m4=2;
M=length(a);
jljl_min=zeros(1,M);
for n=1:M
jljl(:)=sqrt((jlzx(1,:)-maizu(1,n)).^2+(jlzx(2,:)-maizu(2,n)).^2);
jljl_min(n)=min(jljl);
if jljl_min(n)＞500
continue;
end
if jljl(1)==jljl_min(n)
jlzl(1,m1)=n;
m1=m1+1;
elseif jljl(2)==jljl_min(n)
jlzl(2,m2)=n;
m2=m2+1;
elseif jljl(3)==jljl_min(n)
jlzl(3,m3)=n;
m3=m3+1;
else
jlzl(4,m4)=n;
m4=m4+1;
end
end
Mjl(1)=m1-1;Mjl(2)=m2-1;Mjl(3)=m3-1;Mjl(4)=m4-1; %聚类种类每类个数
%对每类求出新的聚类中心
jlzx_n=zeros(2,4);
for n=1:4
sum1=0;sum2=0;
for nn=1:Mjl(n)
sum1=sum1+maizu(1,jlzl(n,nn));
sum2=sum2+maizu(2,jlzl(n,nn));
jlzx_n(1,n)=sum1/Mjl(n);
jlzx_n(2,n)=sum2/Mjl(n);
end
end
%判决条件（新旧两个聚类中心）距离较小
panjue_js=sqrt(sum(sum((jlzx-jlzx_n).^2)));
%设置新的判决中心
jlzx=jlzx_n;
end