search_wide_signal_all.m

蚁群算法的程序

% 基本蚁群算法应用于搜索信号有效带宽
% 信噪比SNR=-10dB
% 信号有效带宽1000~1200Hz
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 第一部分 
% 高斯噪声环境下宽带矢量信号
% 固定信号源，水听器位于原点
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
clc;

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

%信号源方位
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);

%产生噪声
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);

l_length=length(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);
%     I(i)=sqrt(enegry(i,1).^2);
end
%计算路径的长度
l=abs(10*log10(2000*(1./I)));

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

 %总搜索次数为time
time=30;           
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 第四部分
%数据的统计和作图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g=round(l_length/(freq_down-freq_up+1));
%把被检测的信号频率范围划分为g个频带。
%每个频带大小为200Hz 
%统计每个频带上的蚂蚁数量,存于矩阵n_1
for c=1:time
    for d=1:g
        n_1(d,c)=sum(n(100*2*(d-1)+1:100*2*d,c));
    end
end

%找出蚂蚁数量最多的子频带,输出它的序号d_3
for d_1=1:g
    sum_1(d_1)=sum(n_1(d_1,:))/time;
end
Max=sum_1(1);
for d_2=2:g
    if Max<=sum_1(d_2)
       Max=sum_1(d_2); 
       d_3=d_2; 
   end
end

%在循环搜索后统计被检测信号的每个频率平均蚂蚁数目
%输出检测到的信号的有效频带
A(number)=0;
for d_4=1:number
    sum_2(d_4)=sum(n(d_4,:));
    average(d_4)=sum_2(d_4)/time;
    if average(d_4)>8
        A(d_4)=d_4;
    end 
end
for d_5=1:number
    if A(d_5)~=0
        freq_up_1=A(d_5);
        break
    end
end
for d_6=number:-1:1
    if A(d_6)~=0
        freq_down_1=A(d_6);
        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);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
plot(f,10*log10(En_p),'k');
% title('声压功率谱');
xlabel('频率(Hz)');
ylabel('能量(dB)');

figure(2);
plot(f,l,'k');
% title('路径');
xlabel('频率(Hz)');
ylabel('长度(cm)');

figure(3);
plot(n_1(d_3,:),'k');
% title('有效路径(带宽)上的蚂蚁数目');
xlabel('选择次数(次)');
ylabel('蚂蚁数目(只)');

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

n
n_1
d_3
freq_up_1
freq_down_1
SNR0
toc;