一、设计要求
某单波束测深仪最大测量水深为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/sf0 = 15e3; %信号频率 Unit：Hzfs = f0*10; %采样率 最高频率的5倍 Unit：HzTao = 5/f0; %信号脉宽 Unit：sSample_time = 0.1; %采样时长 Unit：s 假设信号发射的时刻为零时刻SNR = 40; %信噪比 Unit：dBH = 10; %水深 Unit：mH1 = 5; %发射换能器水深 Unit：mH2 = 5; %接收换能器水深 Unit：mD = 1; %接收与发射换能器水平距离 Unit：mRe_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);运行后得到图像：
得到理想情况下的仿真数据。
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;enddelay_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 pathsfor 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);endambientnoise_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 PSDnpsd_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 signaltxsig_flip=conj(tx_source(end:-1:1));%matched-filtering%division by ttl_samples is necessary for normalizationest_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 IDenv_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 dBmaxamp_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效果实现如下：
2.3绘制冲激响应图
plotting.m
clearclcfilename = '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 onxlabel('相对时延/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 onxlabel('相对时延/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 onendhold offgrid oncolorbarxlabel('时延(sec)')ylabel('Range(m)')title(filename)生成结果如下：
三.滤波器设计方案
对于本次问题采样IIR滤波器：
切比雪夫Ⅱ型滤波器，种类为带通滤波器
Fstop1=180KHz
Fpass1=190KHz
Fpass2=210KHz
Fstop2=220KHz
Astop1=90
Apass=1
Astop2=90
阶：22阶
节：11
稳定滤波幅值响应
FS：采样频率为2000KHz
滤波器的幅值响应图像：
滤波器设计代码如下：
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 FrequencyFstop1 = 180; % First Stopband FrequencyFpass1 = 190; % First Passband FrequencyFpass2 = 210; % Second Passband FrequencyFstop2 = 220; % Second Stopband FrequencyAstop1 = 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]
滤波后结果如图：
学习不易，转载请联系博主。