【Matlab】海底声学模拟（Bellhop）以及滤波器的设计-Toy模板网

这篇具有很好参考价值的文章主要介绍了【Matlab】海底声学模拟（Bellhop）以及滤波器的设计。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

一、设计要求

某单波束测深仪最大测量水深为300米，请根据《水声学原理》和《数字信号处理》相关知识，仿真设计该单波束测深仪的数字信号处理系统（包括模拟滤波器参数、采样频率、量化精度等工作参数；FIR/IIR滤波器设计，并对数字信号进行：匹配滤波；底检测；底跟踪和声呐图绘制等处理）。

（PS：需要全部代码文件文件请点击这里，需要Bellhop使用说明书请点击这里。）

二、采样数据模拟生成

1.理想条件下声呐采样波形生成

1.1假设出的理想条件：

（1）基于射线声学理论

（2）几何衰减按球面波传播衰减规律衰减，不考虑吸收衰减

（3）仅考虑水底的反射

（4）考虑在高斯白噪声背景下

（5）整个空间声速分布均匀

1.2在假设信号发射的时刻为零时刻的前提下，输入的参数说明

表1 输入参数说明

C=1500	声速	单位：m/s
f0=15e3	信号频率	单位：Hz
fs=f0*10	采样率，最高频率的5倍	单位：Hz
Tao=5/f0	信号脉宽	单位：s
Sample_time=0.1	采样时长	单位：s
SNR=60	信噪比	单位：dB
H=10	水深	单位：m
H1=5	发射换能器水深	单位：m
H2=5	接收换能器水深	单位：m
D=1	接收与发射换能器水平距离	单位：m
Re_coef_bottom	水底反射系数	单位：无

1.3生成发射信号并模拟反射环境

data_generation.m

function receive_signal = data_generation(c,fs,f0,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_bottom)
%**********************************通信声呐水听器数据生成函数************************************
%几点假设：
%（1）基于射线声学理论
%（2）几何衰减按球面波传播衰减规律衰减，不考虑吸收衰减
%（3）仅考虑水底的反射
%（4）考虑在高斯白噪声背景下
%（5）整个空间声速分布均匀
%***************************************输入参数说明******************************************
% c                  声速             Unit：m/s
% fs                 采样率           Unit：Hz
% f0                 信号频率         Unit：Hz
% Tao                信号脉宽         Unit：s
% Sample_time        采样时长         Unit：s      假设信号发射的时刻为零时刻
% SNR                信噪比           Unit：dB
% H                  水深                    Unit：m
% H1                 发射换能器水深           Unit：m
% H2                 接收换能器水深           Unit：m
% D                  接收与发射换能器水平距离  Unit：m
% Re_coef_bottom     水底反射系数  存在半波损失
Ts = 1/fs;      %采样时间间隔
sample_num = fix(Sample_time*fs);       %采样总点数
nTs = (0:sample_num-1)/fs;              %离散的采样时刻
sample_num1 = fix(Tao*fs);              %信号的采样点数
nTs1 = (0:sample_num1-1)/fs;            %信号的离散的采样时刻
receive_signal = zeros(size(nTs));      %用于存储水听器接收的信号
d0 = sqrt((H1-H2)^2+D^2);     %两个换能器的直线距离
S0 = 1/d0;                    %直达波的声压幅值
Noise_var = 10^(-SNR/20)*S0;  %白噪声的方差
%% 计算到达回波信号
real_time = d0/c;                                     %信号的实际到达时刻
signal_start_time = fix(real_time*fs)+1;              %信号第一个采样点的时刻
phase = (signal_start_time*Ts-real_time)/Ts*2*pi;     %得到信号的第一个采样点的相位 
receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) = S0*sin(2*pi*f0*nTs1+phase); %模拟采样
receive_signal = receive_signal + Noise_var*randn(size(nTs));     %加入噪声
M=power(D/(2*H-H1-H2),2);
d1=sqrt(M+1)*(H-H1);
d2=sqrt(M+1)*(H-H2);
dib=  d1+d2  %反射波总路程
Sib=1/dib*Re_coef_bottom;   %反射波的声压幅值
real_time = dib/c;                                                   %信号的实际到达时刻
signal_start_time = fix(real_time*fs)+1;                             %信号第一个采样点的时刻
phase = (signal_start_time*Ts-real_time)/Ts*2*pi;                    %得到信号的第一个采样点的相位 
receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) =...
receive_signal(signal_start_time:(signal_start_time+sample_num1-1)) + Sib*sin(2*pi*f0*nTs1 + phase); %模拟采样                  
receive_signal = receive_signal/S0;    %将直达波的幅值归一化

调用上述函数，主运行函数（test.m）

close all;
clear all;
clc;
c = 1500;            %声速             Unit：m/s
f0 = 15e3;           %信号频率         Unit：Hz
fs = f0*10;          %采样率 最高频率的5倍   Unit：Hz
Tao = 5/f0;         %信号脉宽         Unit：s
Sample_time = 0.1;   %采样时长         Unit：s      假设信号发射的时刻为零时刻
SNR = 40;            %信噪比           Unit：dB
H = 10;       %水深                    Unit：m
H1 = 5;       %发射换能器水深           Unit：m
H2 = 5;       %接收换能器水深           Unit：m
D = 1;       %接收与发射换能器水平距离  Unit：m
Re_coef_bottom = 1;   %水底反射系数  
%%
sample_num = fix(Sample_time*fs);       %采样总点数
nTs = (0:sample_num-1)/fs;              %离散的采样时刻
%生成水听器接收数据
receive_signal0 =data_generation(c,fs,f0,Tao,Sample_time,SNR,H,H1,H2,D,Re_coef_bottom);
figure;
plot(nTs,receive_signal0);
xlabel('时间s');
ylabel('幅度');
string = ['水听器接收的仿真数据,信噪比SNR=',num2str(SNR),'dB'];
title(string);

运行后得到图像：

bellhop模型,海底模拟仿真实验,人工智能

得到理想情况下的仿真数据。

2.Bellhop射线模型与参数设计

2.1Bellhop模型的介绍：

BELLHOP 模型是通过高斯波束跟踪方法(Porter 和 Bueker，1987 年)，计算水平非均匀环境中的声场。高斯波束跟踪方法对于高频水平变化问题特别有吸引力，这是简正波、波数积分和抛物线模型不可替代的。高斯束射线跟踪法的基本思想是将高斯强度分布与每条声线联系起来，该声线为高斯声束的中心声线, 这些声线能较平滑的过渡到声影区，也能较平滑的穿过焦散线，所提供的结果与全波动模型的结果更为一致。

2.2Bellhop模型的使用：

（1）确定env文件。

具体操作请参考Bellhop使用说明书。

（2）设计程序

设计程序实现任意信号输入，实现接收信号的Bellhop模型模拟。

bharr2ir.m

function ir_vec=bharr2ir(delay, amp, sampling_rate)

valid_delay_index=find(delay>0);
if isempty(valid_delay_index)
    disp('[bharr2ir]Error: Zero path simulated by BELLHOP.');
    return;
end
delay_vec=delay(valid_delay_index);
amp_vec=amp(valid_delay_index);

%单位冲激响应进行采样
delay_min=min(delay_vec);
delay_max=max(delay_vec);

cir_length=round(delay_max*sampling_rate);
ir_vec=zeros(cir_length, 1);

%find individual ray paths
for icn=1: length(delay_vec)
    %calculate the arrival index
    %init_delay gives some zeros priror to the first path
    arr_id=round(delay_vec(icn)*sampling_rate);
    
    %generate impulse response. Note that sometime, multiple returns can be
    %generate for the same delay in BELLHOP. 
    ir_vec(arr_id)=ir_vec(arr_id)+amp_vec(icn);
end

ambientnoise_psd.m

function [npsd_db]=ambientnoise_psd(windspeed, fc)
%考虑海洋湍流和风力

%Turbulance noise湍流噪声
ANturb_dB=17-30*log10(fc/1000);

%Ambient noise in dB (wind driven noise)
ANwind_dB=50+7.5*sqrt(windspeed)+20*log10(fc/1000)-40*log10(fc/1000+0.4);

%Thermo noise in dB (wind driven noise)热噪声
ANthermo_dB=-15+20*log10(fc/1000);

%Total noise PSD
npsd_db=10*log10(10^(ANturb_dB/10)+10^(ANwind_dB/10)+10^(ANthermo_dB/10));

signal_mf.m

function est_ir_vec=signal_mf(rx_signal, tx_source, cir_length)

ttl_samples=length(tx_source);

%generate matched-filter from the source signal
txsig_flip=conj(tx_source(end:-1:1));

%matched-filtering
%division by ttl_samples is necessary for normalization
est_ir_vec=conv(txsig_flip, rx_signal)/ttl_samples;

uw_isi.m

function [y, adj_ir_vec]=uw_isi(ir_vec, sl_db, tx_source, nv_db)
%与单位冲激响应函数卷积并加入噪声
%发送信号与信道冲激响应进行卷积
tx_sig=conv(tx_source, ir_vec);

%接收信号信噪比
SNR=sl_db-nv_db;

%信道输出
y=awgn(tx_sig,SNR);
adj_ir_vec=ir_vec;

Bellohop.m

function[rx_signal,Ts]=Bellohop()%rx_signal信号序列 fs采样频率
clear; close all; clc; 

%采样率
sampling_rate=2000e3;

%发射机中心频率
%该值和计算噪声有关
fc=200e3;
windspeed=5;

radio=sampling_rate/fc  %采样率与中心频率的比值
%发射机参数
sl_db=200;  %发射机声源级
bw=8e3;     %带宽

%通过Actup的arr文件获取所需的幅值和时延（单位冲激响应）
%BELLHOP run ID
env_id='env';
%read BELLHOP arr file:
[ amp1, delay1, SrcAngle, RcvrAngle, NumTopBnc, NumBotBnc, narrmat, Pos ] ...
    = read_arrivals_asc( [env_id '.arr'] ) ;
[m,n]=size(amp1);
amp=amp1(m,:);
delay=delay1(m,:);

%Step 1: 创建任意波形（这里以正弦波为例）
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F_Serial_Signal=[zeros(1,1),ones(1,1),zeros(1,10)];  %待发送串行数据
Signal_L=length(F_Serial_Signal);                    %发送数据长度
communication_rate=40e3;                       %通信速率 
Communication_radio= sampling_rate/communication_rate; %采样倍数
signal=repmat(F_Serial_Signal,Communication_radio,1);
signal2=reshape(signal,1,Signal_L*Communication_radio);      %调整后的数据
signal_length=length(signal2);                 %数据长度
t=0:1/sampling_rate:(signal_length-1)/sampling_rate;      %时间
modulation_signal=cos(2*pi()*fc*t); %载波信号
          


tx_source=signal2.*modulation_signal;          %发送

%Step 2:加入噪声和多径影响
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Narrmx=10; %limit ourselves to use the first Narrmax paths

%对单位冲激响应进行采样
ir_vec=bharr2ir(delay, amp, sampling_rate);
%strongest tap in dB
maxamp_db=20*log10(max(abs(amp)))+sl_db;%声源级（dB）加由于信道衰减的能量（dB），表示接收到的信号的能量

%噪声级
[npsd_db]=ambientnoise_psd(windspeed, fc);%npsd_db是噪声的能量谱密度（单位频带内的能量）转换成dB，
nv_db=npsd_db+10*log10(bw);%相当于能量谱密度乘以带宽，即该频带内的噪声的能量
% maxamp_db和nv_db的差值应该为信噪比
disp(['Strongest tap strength=' num2str(maxamp_db,'%.1f') ' dB; Noise variance=' num2str(nv_db,'%.1f') 'dB']);

%考虑多径和噪声影响后的接收信号响应
[rx_signal, adj_ir_vec]=uw_isi(ir_vec, maxamp_db, tx_source, nv_db);

Ts=1/sampling_rate;

t2=0:Ts:(length(rx_signal)-1)*Ts;

figure(1)
subplot(2,1,1)
plot(t,tx_source)
title('发送信号')
xlabel('时间/s')
ylabel('幅值(V)')
set(gca, 'Fontsize', 16);

subplot(2,1,2)
plot(t2,rx_signal)

title('接收到的信号')
xlabel('时间/s')
ylabel('幅值(V)')
set(gca, 'Fontsize', 16);
end

效果实现如下：

bellhop模型,海底模拟仿真实验,人工智能

2.3绘制冲激响应图

plotting.m

clear
clc
filename = 'env.arr';
Minimum_range=100  %（接收水听器的水平方向上接收范围最小值，m）----R 
Maximum_range=1000 %（接收水听器的水平方向上接收范围最大值，m）---RB 

[ amp1, delay, SrcAngle, RcvrAngle, NumTopBnc, NumBotBnc, narrmat, Pos ]... 
 = read_arrivals_asc(  filename );
%%单位冲激响应
[m,n]=size(amp1);
amp = abs(amp1); %取模  
x = delay(m,:); %获取第50个接收机的时延和幅值
y = amp(m,:);
figure(1)
stem(x,y)
grid on
xlabel('相对时延/s')
ylabel('幅度')
title('单位冲激响应')

%%归一化冲激响应
Amp_Delay = [x;y];
Amp_Delay(:,all(Amp_Delay==0,1))=[]; %去掉0值
Amp_Delay=sortrows(Amp_Delay',1);  %按照时延从小到大排序
normDelay = Amp_Delay(:,1)-Amp_Delay(1,1);%归一化时延
normAmp = Amp_Delay(:,2)/Amp_Delay(1,2);%归一化幅度
figure(2)
stem(normDelay,normAmp,'^')
grid on
xlabel('相对时延/s')
ylabel('归一化幅度')
title('归一化冲激响应')

figure(3)
mum=1:m;
ReRange = Minimum_range+(Maximum_range-Minimum_range)/m*mum;
for i=1:min(narrmat)
plot(delay(:,i),ReRange,'o')
hold on
end
hold off
grid on
colorbar
xlabel('时延(sec)')
ylabel('Range(m)')
title(filename)

生成结果如下：

bellhop模型,海底模拟仿真实验,人工智能

三.滤波器设计方案

对于本次问题采样IIR滤波器：

切比雪夫Ⅱ型滤波器，种类为带通滤波器

Fstop1=180KHz

Fpass1=190KHz

Fpass2=210KHz

Fstop2=220KHz

Astop1=90

Apass=1

Astop2=90

阶：22阶

节：11

稳定滤波幅值响应

FS：采样频率为2000KHz

滤波器的幅值响应图像：

bellhop模型,海底模拟仿真实验,人工智能

滤波器设计代码如下：

function Hd = filter_2nd
%FILTER_2ND 返回离散时间滤波器对象。
% MATLAB Code
% Generated by MATLAB(R) 9.11 and DSP System Toolbox 9.13.
% Generated on: 14-Dec-2022 16:14:37
% Chebyshev Type II Bandpass filter designed using FDESIGN.BANDPASS.
% All frequency values are in kHz.
Fs = 2000;  % Sampling Frequency
Fstop1 = 180;         % First Stopband Frequency
Fpass1 = 190;         % First Passband Frequency
Fpass2 = 210;         % Second Passband Frequency
Fstop2 = 220;         % Second Stopband Frequency
Astop1 = 90;          % First Stopband Attenuation (dB)
Apass  = 1;           % Passband Ripple (dB)
Astop2 = 90;          % Second Stopband Attenuation (dB)
match  = 'stopband';  % Band to match exactly
% Construct an FDESIGN object and call its CHEBY2 method.
h  = fdesign.bandpass(Fstop1, Fpass1, Fpass2, Fstop2, Astop1, Apass, ...
                      Astop2, Fs);
Hd = design(h, 'cheby2', 'MatchExactly', match);

% [EOF]

滤波后结果如图：

bellhop模型,海底模拟仿真实验,人工智能