海杂波分析的matlab程序

function [WB SumOfAll fd FirstOrder]=AllOrder(f0,Vp,Q,QN,Q0,Reso) %海浪散射系数仿真程序 个参数下面有说明 可以计算单基地雷达和双基地雷达的海浪散射系数 DopVal=zeros(1,401); % f0=15; Vp=15; Q=0*pi/180;QN=90*pi/180.0;Q0=0*pi/180.0;Reso=2000; g=9.8; A=4/(3*pi); s=4; K0=2*pi*f0/300; %//无线电波长 WB=sqrt(2*g*K0*cos(Q0)); % //Bragg 基频 %/* 参数声明 % Vp---输入风速19.5(m/s) % f0---雷达操作频率 % K0---雷达波数 % K1---一阶海洋波数 % K2---二阶海洋波数 % Wd---多普勒频率 % WB---Bragg多普勒频率 % Q ---风向（角度) % QN---散射椭圆法线矢量方向 % Q0---半双基地角 % Q1---一阶海洋波方向 % Q2---二阶海洋波方向 % s ---扩展因子 */ % //计算二阶散射截面积 % //-3WB<Wd<WB部分，步长0.025WB Singular1=0; Singular2=0; Singular3=0; Singular4=0; Singular5=0; Singular6=0; Singular7=0; Singular8=0; Singular9=0; Singular10=0; for N=1:160 Wd=(-5.0+0.025*(N-1))*WB; if(Wd+sqrt(2.0)*WB>0) %// Wd>=-sqrt(2)*WB部分 if(N==145) Wd=-sqrt(2)*WB; %// 奇异值 Singular1=N; end if(abs(Wd+(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB) Wd=-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB; Singular2=N; %//电场耦合系数产生的奇异值 end DopVal(N)=dif(QN-pi,QN,1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN,QN+pi,1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0); else %// Wd<-sqrt(2)WB if(abs(Wd+(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB) Wd=-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB; Singular3=N; %//电场耦合系数产生的奇异值 end if(abs(Wd+(2^0.75)*sqrt(((1+sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB) Wd=-(2^0.75)*sqrt(((1+sin(Q0))^0.5)/cos(Q0))*WB; Singular4=N; end DopVal(N)=dif(QN-pi,QN-acos(2*(WB/Wd)*(WB/Wd)),1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN+acos(2*(WB/Wd)*(WB/Wd)),QN+pi,1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0); end end % //计算-WB<Wd<0部分 for N=162:200 Wd=(-5+0.025*(N-1))*WB; if(abs(Wd+(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB) Wd=-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB; Singular5=N; end DopVal(N)=dif(QN-pi,QN,-1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN,QN+pi,-1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0); end N=N+1; if(N==201) % //Wd=0 Wd=(-5+0.025*(N-1))*WB; DopVal(N)=dif(QN-0.5*pi,QN+0.5*pi,-1,1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0); end for N=202:240 Wd=(-5+0.025*(N-1))*WB; if(abs(Wd-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB) Wd=(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB; Singular6=N; %//电场耦合系数产生的奇异值 end DopVal(N)=dif(QN-pi,QN,1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN,QN+pi,1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0); %// printf("%d,%e\n",N,DopVal[N]); end %//计算WB<Wd<3WB部分 for N=242:401 Wd=(-5+0.025*(N-1))*WB; if(Wd-sqrt(2.0)*WB<0) %// Wd>=sqrt(2)*WB部分 if(N==257) Wd=sqrt(2)*WB; %// 奇异值 Singular7=N; end if(abs(Wd-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB) Wd=(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB; Singular8=N; end %//电场耦合系数产生的奇异值 DopVal(1,N)=dif(QN-pi,QN,-1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN,QN+pi,-1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0); else % //Wd>sqrt(2)*WB if(abs(Wd-(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB) Wd=(2^0.75)*sqrt(((1-sin(Q0))^0.5)/cos(Q0))*WB; Singular9=N; end if(abs(Wd-(2^0.75)*sqrt(((1+sin(Q0))^0.5)/cos(Q0))*WB)<=0.0125*WB) Wd=(2^0.75)*sqrt(((1+sin(Q0))^0.5)/cos(Q0))*WB; Singular10=N; end DopVal(N)=dif(QN-pi,QN-acos(2*(WB/Wd)*(WB/Wd)),-1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0)+dif(QN+acos(2*(WB/Wd)*(WB/Wd)),QN+pi,-1,-1,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0); end end DopVal(201)=(DopVal(200)+DopVal(202))/2; SecondDopVal=DopVal; SecondDopVal1=DopVal; SecondDopVal=SmoothScecondOrderData(SecondDopVal,Singular1,Singular2,Singular3,Singular4,Singular5,Singular6,Singular7,Singular8,Singular9,Singular10); %计算一阶海杂波 FirstOrder=zeros(1,401); for N=1:401 if(N-1>200) m=-1; else m=1; end K=((-5+0.025*(N-1))*WB)^2/g; if(N-1==160|N-1==240) Val=(2^4)*pi*(K0^2)*0.00405*A*(K^-4)*exp(-0.74*g*g*(K^-2)*(Vp^-4))*(cos(((1-m)/2*pi+QN-Q)/2.0)^s)*(K^2.5)*cos(Q0)/sqrt(g)*Reso*(cos(Reso/2*(K/cos(Q0)-2*K0))^2); else Val=(2^4)*pi*(K0^2)*0.00405*A*(K^-4)*exp(-0.74*g*g*(K^-2)*(Vp^-4))*(cos(((1-m)/2*pi+QN-Q)/2.0)^s)*(K^2.5)*cos(Q0)/sqrt(g)*Reso*(sin(Reso/2*(K/cos(Q0)-2*K0))/(Reso/2*(K/cos(Q0)-2*K0))^2); end if(Val>=FirstOrder(N)) FirstOrder(N)=Val; else end end FirstOrder(201)=(FirstOrder(200)+FirstOrder(202))/2; FirstOrder=dBConverter(FirstOrder); FirstOrder=SmoothData(FirstOrder); FirstOrder=10.^(FirstOrder/10); %给出一阶和二阶海杂波RCS的和 SumOfAll=SecondDopVal+FirstOrder; %给出多普勒坐标fd N=1:401; Wd=(-5.0+0.025*(N-1))*WB; fd=Wd; %plot(fd,10*log10(SumOfAll)); %%定义求数组转化成dB表示的函数 function dBVal=dBConverter(AA) for L=1:length(AA) if(AA(L)<=0.000000) AA(L)=-90.0; else AA(L)=10*log10(AA(L)); end if(AA(L)<-90.0) AA(L)=-90.0; end end dBVal=AA; %定义平滑曲线的函数 function SmoothVal=SmoothData(B) %求极大值点，比把极大值点坐标保存在jj数组里面，之后把极大值点数据保存在flav里面 jj=[]; for j=2:length(B)-1 if((B(j)>(B(j-1)))&(B(j)>B(j+1))) jj(end+1)=j; end end jj=[1 jj length(B)]; jjj=length(jj); flatv=[]; for i=1:jjj flatv(i)=B(jj(i)); end %插值进行曲线平滑 i=1:length(B); SmoothVal=interp1(jj,flatv,i,'spline'); %定义功能函数平滑二阶数据 function OutSmoothResult=SmoothScecondOrderData(InputSecondOderData,Singular1,Singular2,Singular3,Singular4,Singular5,Singular6,Singular7,Singular8,Singular9,Singular10) NN=length(InputSecondOderData); XX=InputSecondOderData; for i=1:NN if(XX(i)<=0) XX(i)=-90; else XX(i)=10*log10(XX(i)); end if(XX(i)<-90) XX(i)=-90; end end for i=2:NN-1 if(i==Singular1|i==Singular2|i==Singular3|i==Singular4|i==Singular5|i==Singular6|i==Singular7|i==Singular8|i==Singular9|i==Singular10) if(XX(i)>XX(i-1)+15) XX(i)=XX(i-1)+15; end else if((XX(i)>XX(i-1)&XX(i)>XX(i+1))&(XX(i)>(max(XX(i-1),XX(i+1))+2))) XX(i)=(XX(i-1)+XX(i+1))/2; end end end XX=10.^(XX./10); OutSmoothResult=XX; %定义求二阶海杂波积分的函数 function val=dif(LowLimit,UpLimit,m1,m2,Wd,WB,f0,Vp,Q,QN,Q0,Reso,K0) C=3.0e8; %光速 g=9.8; A=4/(3*pi); s=4; val=0; Sigma2=0.0; %初始化积分值 n3=m1*m2; if(abs(Wd)>sqrt(2)*WB) Beta=acos(2*(WB/Wd)*(WB/Wd)); else Beta=0; end h=0.5*pi/180; %Theta=0.5的积分步长 Q1=LowLimit; %积分下限 while(Q1<=UpLimit) Sigma=0.0; %初始积分值 % 在给定Wd,Q1时，用牛顿法计算海洋波数K1*/ K_sqrt=0.00000001; %猜测的K1的初始均方根值 fK_sqrt=10; while((abs(fK_sqrt))>=0.000001) if(K_sqrt<0) K_sqrt=-K_sqrt; else K_sqrt=K_sqrt; end fK_sqrt=Wd+sqrt(g)*(m1*K_sqrt+m2*((((K_sqrt^4)+((2*K0*cos(Q0))^2)-4*K_sqrt*K_sqrt*K0*cos(Q0)*cos(Q1-QN)))^0.25)); % fK_sqrt=Wd+m1*sqrt(g*K1)+m2*sqrt(g*K2) fK_sqrtP=sqrt(g)*(m1+m2*K_sqrt*(K_sqrt*K_sqrt-2*K0*cos(Q0)*cos(Q1-QN))*((((K_sqrt^4)+((2*K0*cos(Q0))^2)-4*K_sqrt*K_sqrt*K0*cos(Q0)*cos(Q1-QN)))^-0.75)); % fK_sqrtP 是fK_sqrt的微分 K_sqrt=K_sqrt-fK_sqrt/fK_sqrtP; %K_sqrt的新值 end K1=K_sqrt*K_sqrt; %一阶海洋波