search_wide_signal_1.m

蚁群算法的程序

%基本蚁群算法应用于搜索信号有效带宽
%信噪比SNR=-10dB
%信号有效带宽250~300Hz

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%高斯噪声环境下宽带矢量信号
%固定信号源，水听器位于原点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
clc;

%初始化参数
x_a=0;
y_a=0;           %水听器坐标
len=8192;        %采样点数
fs=4096;         %采样频率
SNR=-10;         %信噪比
nfft=4096;
%信号带宽
freq_up=301;     
freq_down=350;
%待检测频率范围
up=201;
down=500;

%信号源方位
target_x=30;
target_y=40;
r_s=sqrt((target_x-x_a)^2+(target_y-y_a)^2); 
phi_s=atan2(target_y-y_a,target_x-x_a);
ini_phi_s=phi_s*180/pi;

%产生噪声
An=8;                  %噪声幅度基值
W=randn(len,4);
N=[1,0,0,0;
   0,1/3,0,0;
   0,0,1/3,0;
   0,0,0,1/3];
N=sqrt(N);
noise=W*N;
p_n=noise(:,1);
En_n=std(An*p_n).^2;            %噪声能量	
df=fs/2;           
dEn_n=En_n/df;                  %能量密度



%产生随机信号
signal_p=randn(len,1);
%信号通过椭圆带通滤波器
wn=[freq_up freq_down]*2/fs;
[b,a]=ellip(6,0.1,30,wn);
signal=filter(b,a,signal_p);

%带宽信号
En_s=std(signal).^2;                    %信号能量
dEn_s=En_s/(freq_down-freq_up);         %信号能量密度
A_s=sqrt(2*10.^(SNR/10)*En_n/En_s);   %信号幅值
p_s=A_s.*signal;       %声压  
vx_s=p_s.*cos(phi_s);  %X方向振速
vy_s=p_s.*sin(phi_s);  %Y方向振速

%信号+噪声
p=p_s+An*noise(:,1);
En=sum(p.^2)/length(p);
vx=vx_s+An*noise(:,2);
vy=vy_s+An*noise(:,3);
vz=An*noise(:,4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%建立信号频率与蚁群算法中路径长度的关系
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%功率谱估计
window=hanning(nfft);               %加汉宁窗
[En_p,f]=psd(p,nfft,fs,window);
[En_vx,f]=psd(vx,nfft,fs,window);
[En_vy,f]=psd(vy,nfft,fs,window);
[pvx,f]=csd(p,vx,nfft,fs,window);
[pvy,f]=csd(p,vy,nfft,fs,window);
[pvz,f]=csd(p,vz,nfft,fs,window);

figure(1);
plot(f,10*log10(En_p),'k');
title('声压功率谱');
xlabel('frequence/Hz');
ylabel('Energy/dB');

En_p=En_p(up:down);
En_vx=En_vx(up:down);
En_vy=En_vy(up:down);
pvx=pvx(up:down);
pvy=pvy(up:down);
pvz=pvz(up:down);

% %估计方位
% En_r=En_vx+En_vy;
% energy=real(pvx).^2+real(pvy).^2;			%求出各个频率能量
% 
% estimate_theta=atan2(real(pvy),real(pvx))*180/pi;				%求出各个频率指定的方位
% amplitude_theta=energy./(En_p.*En_r).*(En_p+En_r);            %各谱线幅度信息
% ind_max=find(amplitude_theta==max(amplitude_theta));
% 
% result_phi_s=estimate_theta(ind_max);      %能量最大的方位信息

temp_f=f(up:down);
l_length=length(temp_f);  
for i=1:l_length
    enegry(i,1)=real(En_p(i));
    enegry(i,2)=real(pvx(i));
    enegry(i,3)=real(pvy(i));
    enegry(i,4)=real(pvz(i));
    I(i)=sqrt(enegry(i,2).^2+enegry(i,3).^2+enegry(i,4).^2);
end
l=abs(10*log10(1000./I))+2;

figure(2);
plot(f(up:down),l,'k');
title('路径');
xlabel('frequence/Hz');
ylabel('length');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%蚁群算法
%参考文献：《蚁群算法原理及其应用》 段海滨
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%参数初始化
number=down-up+1;
q=1:number;
alpha=0.1;      %信息启发式因子
beta=0.1;       %期望启发式因子
tao(q)=0;        %路径的初始信息量
rou=0.5 ;        %信息素挥发系数
N=2000;         %蚁群大小
Q=30;          %信息素强度
D(q)=l;       % 路径的长度设定
yita=1./l';        %路径的启发函数

  %总搜索次数为time
time=20;           
n(q,time)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%路径上信息素为零，蚁群第一次对路径进行随机选择
for i_0=1:N
    p=rand;
    p_0=1/number;
    sum_p_0=0;
    for j_0=1:number
        sum_p_0=sum_p_0+p_0;
        if sum_p_0>p
            n(j_0,1)=n(j_0,1)+1;
            break
        end
    end
end
            
%第一次选择后，路径上信息素改变
detatao=n(: ,1).*Q./l';
tao=(1-rou).*tao'+detatao;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%蚁群从第二次开始连续进行time-1次循环选择
for k=2:time
    detatao=0;
    for i=1:N
        %每只蚂蚁在选择路径之前，先计算各条路径的选择概率
        sum_1=sum(tao.^alpha.*yita.^beta);
        p(q)=tao.^alpha.*yita.^beta./sum_1;
        %运用轮盘赌选择算法
        p_1=rand;
        sum_p=0;
        for j=1:number
            sum_p=sum_p+p(j);
            if sum_p>p_1;
                n(j,k)=n(j,k)+1;
                break
            end
        end
    end
    
%信息素更新
detatao=n(:,k).*Q./l';
tao=(1-rou)*tao+detatao;
end

%在循环搜索后统计被检测信号的每个频率平均蚂蚁数目
%输出检测到的信号的有效频带
A(number)=0;
for d=1:number
    sum_a(d)=sum(n(d,:));
    average(d)=sum_a(d)/time;
    if average(d)>10
        A(d)=d;
    end 
end
for d_1=1:number
    if A(d_1)~=0
        freq_up_1=A(d_1)+200;
        break
    end
end
for d_2=number:-1:1
    if A(d_2)~=0
        freq_down_1=A(d_2)+200;
        break
    end
end
%噪声通过椭圆带通滤波器
wn0=[freq_up_1 freq_down_1]*2/fs;
[b0,a0]=ellip(6,0.1,30,wn0);
noise0=filter(b0,a0,p_n);
En_n0=std(An*noise0).^2;
En_s1=std(p_s).^2;  % 信号能量
SNR0=10*log10(En_s1/En_n0);

%估计方位
freq_up_2=freq_up_1-200;
freq_down_2=freq_down_1-200;
En_p=En_p(freq_up_2:freq_down_2);
En_vx=En_vx(freq_up_2:freq_down_2);
En_vy=En_vy(freq_up_2:freq_down_2);
pvx=pvx(freq_up_2:freq_down_2);
pvy=pvy(freq_up_2:freq_down_2);
pvz=pvz(freq_up_2:freq_down_2);
En_r=En_vx+En_vy;
energy=real(pvx).^2+real(pvy).^2;			%求出各个频率能量
estimate_theta=atan2(real(pvy),real(pvx))*180/pi;				%求出各个频率指定的方位
result_phi_s=sum(estimate_theta(:))/(freq_down_1-freq_up_1+1);      %信号平均方位信息


%找出蚂蚁数量最多的子频带,输出它的序号c_6
for k_1=1:time
    for c_3=1:6
        n_3(c_3,k_1)=sum(n(50*(c_3-1)+1:50*c_3,k_1));
    end
end

for c_4=1:6
    n_4(c_4)=sum(n_3(c_4,:))/time;
end
B=n_4(1);
for c_5=2:6
    if B<=n_4(c_5)
        B=n_4(c_5);
        c_6=c_5;
    end
end



figure(3);
plot(n_3(c_6,:),'k');
title('有效路径(带宽)上的蚂蚁数目');
xlabel('number');
ylabel('time');

figure(4);
plot(f(up:down),average,'k');
% title('搜索完成后各个频率上的蚂蚁分布');
xlabel('频率(Hz)');
ylabel('蚂蚁数目(只)');

n
n_3
n_4
freq_up_1
freq_down_1
SNR0
c_6
ini_phi_s
result_phi_s