常规波束形成——时域与频域常规窄带波束形成、信噪比计算(附代码)

这篇具有很好参考价值的文章主要介绍了常规波束形成——时域与频域常规窄带波束形成、信噪比计算(附代码)。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

最近学习了一下《最优阵列处理技术》,应老师要求写一个线性均匀水听器阵列的常规波束形成,由于是初学者,写的可能会有点问题,欢迎大家提出修改建议和指导,写这个主要是记录自己的思考,其次是和初学者进行交流提升。

1. 主要原理阐述

要实现波束形成,首先得了解频率波束响应,波束形成与之息息相关。

1.1 频率波数响应

对于一个水听器阵列,我们要想的得到其对外部信号场的响应,则需要对信号场在空域上进行采样,得到一组信号,表示为:

对每个阵元的输出用一个线性时不变滤波器进行处理,该滤波器的冲激响应为,并对所有输出求和,得到阵列的输出,该过程如下图示:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

针对线性阵列,我们用一张图来说明一下:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

设均匀阵元排列间隔为,假如信号平行射入,信号入射角度与端射方向呈角,则可以通过计算知道,相邻的阵元接收到信号的间隔为:

 设第一个阵元接收到的信号为,则第二个阵元接收到的信号为频域波束形成与时域波束形成区别,matlab,开发语言,矩阵,以此类推,在理想情况下,第n个阵元接收到的信号为:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

为了能够还原我们想要的信号,我们可以使用一个延时求和波束形成器,其原理如下:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

 式子表示为:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

也可以写成频域内的简洁形式:

  其中,称作流形矢量,定义为:

右上角的H为共轭转置的意思。

 其中:

为信号的频率,为方向余弦,为介质中的信号传播速度。

为了进一步讲解,我们假设一个基函数,表示如下:

阵列处理器对一个平面波的响应:

 可以表示为:

在频域内,可以写为:

定义:

 为阵列的频率波束响应。

波束方向图的定义是用入射方向表示频率波数响应函数(在代码中我们使用的角度为60°,可以自行修改)。

2. 代码讲解

2.1 基函数

在代码中,我们沿用基函数:

在我们进行频率波数响应推导的时候,细心的朋友应该发现了,信号接收时间间隔的定义即为的定义,为位置表示,由于是在处理一个平面波信号,方向余弦仅需取一个方向,这里为,故根据的定义:

假设我们的信号为,其中,为信号声源的频率,那么我们现在可以表示每个阵元的接收信号:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

打开有:

该部分代码为:

%% 目标信号生成
f0=300;               % 声源频率,Hz
s=exp(1j*2*pi*f0.*t); % 目标复数信号,1 x LEN  "1j为复数,1为数字一,j为负数"  
as=exp(1j*2*pi*ra*f0*cos(thetas_rad)/c); % 目标驾驶向量,N x 1,驾驶向量=exp(-2*pi*f0*τ*1j*)
snr=40;               % 信噪比设置
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的水听器输出时间波形实数矩阵,Add white Gaussian noise to signal
x = rx.';          % 输出波形矩阵转置,便于二维数组FFT列运算,LEN x N,即接收到的信号

s为源信号,as为流形矢量(也有称驾驶向量的说法),x为得到的输出信号。

注意:正常来讲,我们在进行信号接收时,应该在一定频率内进行扫描,获取不同频率信号再进行筛选;在代码中,我们假设的信号实际是简化了,相当于我们本来就知道信号的频率(设为300Hz),现在只需要在该频率下进行角度扫描,确定信号的来源

2.2 延时处理

接着,我们刚刚提到了延时求和波束形成器,根据频率响应函数的定义,我们需要对信号基函数进行延时求和,由此得到波束方向图:

 我们在1.1频率波数响应函数末尾提到,在代码中我们已经确定了搜索的信号频率,因此仅需要对角度进行扫描,综上所述,代码如下:

%% 目标方位角设置
thetas_deg=60;                 % 目标方位角,度deg
thetas_rad=thetas_deg*deg2rad; % 目标方位角,弧度rad
theta = 0:0.1:180;                  % 任意信号入射角度


%% 频域常规波束形成
tic
for i = 1:length(theta)
    tao = cosd(theta(i))*ra/c;   %在角度theta(i)下得到的时延
    pt_sig = exp(1j*2*pi*f0*tao)' * x.';  %波束,1*LEN
    s_sig(i,:) = pt_sig/N;    %size(s_sig) = i*LEN
end
toc

我们的角度扫描从0到180°,步长为0.1,因此可以得到1801长度的theta,tao为时延补偿,对应中各个阵元的延时补偿。

注意

1. 频率波束响应公式中的转置为共轭转置,在MATLAB中使用单引号,x为输出信号,根据公式不需要转置,这里进行转置主要是为了进行矩阵运算。

2. pt_sig的计算是将每个阵元的对应时延与输出信号相乘得到修正后的信号,然后再相加得到频率响应(频率响应的定义为响应和基函数的卷积累加):

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

如果还是不太明白可以跑出代码后对每个矩阵参数的行列数进行细致观察,从而加深理解。

结果:得到了s_sig,其为i*LEN大小的数据,i为theta的长度。矩阵一列的意义为:每个角度下的频率波数响应;矩阵一行的意义为:同一角度下不同采样点的频率波数响应。

3. 结果展示

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

分析:出现了两个波峰,一个对应60°,一个对应120°,暂且不清楚原因,后续可能会跟进修正,恳请知道原因的朋友赐教,以上。

原因解释:

由于代码基函数部分得到的信号只取了实数部分,后续与修正时延tao的共轭转置相乘时无法通过正交消除120°干扰,如下:

snr=40;               % 信噪比设置
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的水听器输出时间波形实数矩阵,Add white Gaussian noise to signal
x = rx.';          % 输出波形矩阵转置,便于二维数组FFT列运算,LEN x N,即接收到的信号

修改代码:

%% 频域常规波束形成
tic
xx = (as*s).';
for i = 1:length(theta)
    tao = cosd(theta(i))*ra/c;   %在角度theta(i)下得到的时延
    pt_sig = exp(1j*2*pi*f0*tao)' * xx.';  %波束,1*LEN
    s_sig(i,:) = pt_sig/N;    %size(s_sig) = i*LEN
end
toc
plot(theta,s_sig(:,1),'b.-');title("频域常规波束形成");  

得到图像:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

附修改后的MATLAB代码:

clear
close all
clc

% 度与弧度常数
deg2rad = pi/180; % 度转换成弧度
rad2deg = 180/pi; % 弧度转换成度


%% 水听器阵列参数设置
c = 1500;             % 水体声速,米/秒 m/s
N = 70;               % 阵元数
sqrtN=sqrt(N);        % 驾驶向量归一化常数
T=1;                 % 数据时长,秒sec
FS=24e2;              % 采样频率,赫兹Hz
LEN=T*FS;             % 数据长度,点数
t=(0:LEN-1)./FS;      % 数据时间轴,秒sec,1 x LEN 
fc=300;               % 阵列中心频率,赫兹Hz
lambda = c/fc;        % 波长,米m
lambda2 = lambda / 2; % 半波长,米m   d<(lambda/2)?
d = 0.27;             % 阵元间距,米m
ra=d*(0:N-1).';       % 阵元坐标,N x 1   

%% 目标方位角设置
thetas_deg=90;                 % 目标方位角,度deg
thetas_rad=thetas_deg*deg2rad; % 目标方位角,弧度rad
theta = 0:0.1:180;                  % 任意信号入射角度

%% 目标信号生成
f0=300;               % 声源频率,Hz
s=exp(1j*2*pi*f0.*t); % 目标复数信号,1 x LEN  "1j为复数,1为数字一,j为负数"  
as=exp(1j*2*pi*ra*f0*cos(thetas_rad)/c); % 目标驾驶向量,N x 1,驾驶向量=exp(-2*pi*f0*τ*1j*)
snr=40;               % 信噪比设置
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的水听器输出时间波形实数矩阵,Add white Gaussian noise to signal
x = rx.';                                 % 输出波形矩阵转置,便于二维数组FFT列运算,LEN x N,即接收到的信号

%% 频域常规波束形成
tic
xx = (as*s).';
for i = 1:length(theta)
    tao = cosd(theta(i))*ra/c;   %在角度theta(i)下得到的时延
    pt_sig = exp(1j*2*pi*f0*tao)' * xx.';  %波束,1*LEN
    s_sig(i,:) = pt_sig/N;    %size(s_sig) = i*LEN
end
toc
figure(1);
subplot(2,1,1);
plot(theta,s_sig(:,1),'b.-');title("频域常规波束形成");  
subplot(2,1,2);
sig_db = 20*log10(abs(s_sig(:,1).^2));
plot(theta,sig_db,'r.-');title("频域常规波束形成——分贝图");  

如果本篇文章对你有帮助,解决了你的疑惑或者问题,请不要吝啬你的点赞哦,谢谢!

4. 添加频域窄带波束形成方法

clear
close all
clc

% 度与弧度常数
deg2rad = pi/180; % 度转换成弧度
rad2deg = 180/pi; % 弧度转换成度


%% 阵列参数设置
c = 1500;             % 水体声速,米/秒 m/s
N = 35;               % 阵元数
sqrtN=sqrt(N);        % 驾驶向量归一化常数
T=1;                 % 数据时长,秒sec
FS=24e3;              % 采样频率,赫兹Hz
LEN=T*FS;             % 数据长度,点数
t=(0:LEN-1)./FS;      % 数据时间轴,秒sec,1 x LEN 
fc=300;               % 阵列中心频率,赫兹Hz
lambda = c/fc;        % 波长,米m
lambda2 = lambda / 2; % 半波长,米m
d = 0.27;             % 阵元间距,米m
ra=d*(0:N-1).';       % 阵元坐标,N x 1

%% 方位角设置
thetas_deg=90;                 % 目标方位角,度deg
thetas_rad=thetas_deg*deg2rad; % 目标方位角,弧度rad

%% 目标信号生成
f0=3000;               % 声源频率,Hz
f1=f0;                % 起始处理频率,Hz
B=10;                 % 处理带宽,Hz
s=exp(1j*2*pi*f0.*t); % 目标复数信号,1 x LEN
sf=fft(s.',LEN);      % 目标信号FFT,LEN x 1
ks=2*pi*f0/c;         % 声源波数,1/m
% 指向目标方位的驾驶向量生成
as=exp(-1j*2*pi*ra*f0*cos(thetas_rad)/c); % 目标驾驶向量,N x 1
rf=sf*(as.');         % 指向目标方向的输出复数信号,LEN x N
% 指向目标方位的输出信号生成
snr=40;               % 信噪比设置
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的输出时间波形实数矩阵
x = rx.';                                 % 输出波形矩阵转置,便于二维数组FFT列运算,LEN x N
x_fft = fft(x,LEN);   % 对实信号进行FFT

%% 基础数据
df = FS / LEN;        % 采样分辨率,即某个采样点对应的频率
figure(1);plot((0:LEN-1)*FS / LEN,abs(x_fft(:,1)));  % 测试图像,观察目标信号傅里叶变换:频率为f0 = 300
place = round(f0/df) + 1;    % 获取信号FFT频率对应位置
sig = x_fft(place,:);    % 获取频率为f0信号的FFT结果  1*N   频率有0频
theta_N = 0:0.1:180;  % 角度遍历
ntheta = length(theta_N);


%% 常规频域窄带波束
for i = 1:ntheta
    tao = cos(theta_N(i)*deg2rad)*ra/c;  % 在theta_N(i)角度下的时延  N*1
    A = exp(1j*2*pi*f0*tao);
    B = A' * sig.';    % 常规窄带波束  1*N*N*1
    pt_sig(i) = B;
    power_sig(i) =  B^2;  % 功率
    tt(i) = cosd(theta_N(i));
end

pt_abs = abs(pt_sig);
pt2one = pt_abs ./ max(pt_abs); % 归一化

power_abs = abs(power_sig);
power2one = power_abs ./ max(power_abs);  % 归一化
 
figure(2)
plot(tt,20*log10(pt2one),'r');
% plot(theta_N,(pt2one),'r');
grid on

figure(3)
plot(tt,10*log10(power2one),'b');    % 功率的是10
% plot(theta_N,(power2one),'b');
grid on

 运行结果:

 频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

5. 信噪比提升计算

应评论区要求,增加信噪比提升计算代码。该代码基于时域波束形成,主要的思想是:将添加噪声的信号未添加噪声的信号进行波束形成并计算功率,两者相减得到高斯白噪声的功率。最后使用:

 

得到信噪比。

求得角度遍历下的信噪比最大值对应信号所在的方位。

附代码:

clear
close all
clc
 
% 度与弧度常数
deg2rad = pi/180; % 度转换成弧度
rad2deg = 180/pi; % 弧度转换成度
 
%% 水听器阵列参数设置
c = 1500;             % 水体声速,米/秒 m/s
N = 4;               % 阵元数
sqrtN=sqrt(N);        % 驾驶向量归一化常数
T=1;                 % 数据时长,秒sec
FS=24e2;              % 采样频率,赫兹Hz
LEN=T*FS;             % 数据长度,点数
t=(0:LEN-1)./FS;      % 数据时间轴,秒sec,1 x LEN 
fc=300;               % 阵列中心频率,赫兹Hz
lambda = c/fc;        % 波长,米m
lambda2 = lambda / 2; % 半波长,米m   d<(lambda/2)?
d = 0.27;             % 阵元间距,米m
ra=d*(0:N-1).';       % 阵元坐标,N x 1   
 
%% 目标方位角设置
thetas_deg=90;                 % 目标方位角,度deg
thetas_rad=thetas_deg*deg2rad; % 目标方位角,弧度rad
theta = 0:0.1:180;                  % 任意信号入射角度
 
%% 目标信号生成
f0=300;               % 声源频率,Hz
s=exp(1j*2*pi*f0.*t); % 目标复数信号,1 x LEN  "1j为复数,1为数字一,j为负数"  
as=exp(1j*2*pi*ra*f0*cos(thetas_rad)/c); % 目标驾驶向量,N x 1,驾驶向量=exp(-2*pi*f0*τ*1j*)
snr=-20;               % 信噪比设置
SNR_list = zeros(1,100);
Sig = sum(sum(real(as*s).^2,2)/LEN)/N;  % 信号原始功率

tic
for ii = 1:100  % 进行蒙特卡洛模拟测试
rx=awgn(real(as*s), snr, 'measured');     % 含加性高斯噪声的水听器输出时间波形实数矩阵,Add white Gaussian noise to signal 
%% 常规波束形成——不带噪声数据
% tic
xx = real((as*s).');    % 不带噪声的原始数据
for i = 1:length(theta)
    tao = cosd(theta(i))*ra/c;   %在角度theta(i)下得到的时延
    pt_sig = exp(1j*2*pi*f0*tao)' * xx.';  %波束,1*LEN
    Beamforming(i,:) = pt_sig / N; % i * LEN
    s_sig(i,:) = sum((abs(pt_sig/N)).^2)/LEN;    %size(s_sig) = i*LEN
end
% toc

%% 常规波束形成——带噪声数据
% tic
x = rx.';               % 输出波形矩阵转置,便于二维数组FFT列运算,LEN x N,即接收到的信号
for i = 1:length(theta)
    tao = cosd(theta(i))*ra/c;   %在角度theta(i)下得到的时延
    pt_sig2 = exp(1j*2*pi*f0*tao)' * x.';  %波束,1 * N * N * LEN = 1 * LEN
    s_sig2(i,:) = sum((abs(pt_sig2/N)).^2)/LEN;    % i * 1
end
% toc

N_sig = s_sig2 - s_sig; % 噪声功率
SNR_ = 10*log10(s_sig ./ N_sig); % 通过波束形成后数据求信噪比,SNR = 10*log10(信号功率/噪声功率)
% SNR_ = 10*log10(Sig / N_sig); % 通过原始数据功率求信噪比
% 找最大波束处SNR为增益后的信噪比
SNR_bias = max(SNR_) - snr; % 信噪比的差值
SNR_list(ii) = SNR_bias;
end
toc

% 查看波束图
figure;plot(theta,(abs((Beamforming(:,1)))));title("波束形成图")
Anouncement_ = ['SNR = ', num2str(max(SNR_)),'dB,提升为:', num2str(mean(SNR_list)),'dB'];
disp(Anouncement_)

 运行结果展示:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

结果分析:

 可以看到,当阵元数N = 4时,SNR提升大概为6dB左右;后将N = 40,再进行计算,得到结果:

频域波束形成与时域波束形成区别,matlab,开发语言,矩阵

可以看到SNR有所提升,符合增加阵列孔径的引起的信号分辨率提升规律。文章来源地址https://www.toymoban.com/news/detail-758090.html

到了这里,关于常规波束形成——时域与频域常规窄带波束形成、信噪比计算(附代码)的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 信噪比的定义及计算方法

    英文名称叫做SNR或S/N(SIGNAL-NOISE RATIO),又称为讯噪比。是指一个电子设备或者电子系统中 信号与噪声的比例 。这里面的信号指的是来自设备外部需要通过这台设备进行处理的电子信号,噪声是指经过该设备后产生的原信号中并不存在的无规则的额外信号(或信息),并且该

    2024年02月01日
    浏览(35)
  • Matlab 生成一定信噪比的信号

    信噪比公式 1 : S N R = 10 ∗ l o g 10 P s P n 信噪比公式1:SNR=10*log_{10}frac{P_s}{P_n} 信噪比公式 1 : SNR = 10 ∗ l o g 10 ​ P n ​ P s ​ ​ 其中, P s P_s P s ​ 和 P n P_n P n ​ 分别指 信号的平均功率、噪声的平均功率。 当信号和噪声的长度相等为N时,根据平均功率的公式 P s = E s /

    2024年02月11日
    浏览(46)
  • 音频筑基:信噪比SNR指标

    在分析音频信号中,信噪比是我们经常遇到的概念,这里谈谈自己的理解。 定义 SNR,Signal to Noise Ratio,信噪比,也常缩写为S/N 概念 顾名思义,就是信号和噪声的比值,实际应用时比值结果常转到dB域中 信号,原系统输出的理论信号 噪声,指经过某系统后,原信号不存在的

    2024年02月04日
    浏览(42)
  • 视频的动态范围、信噪比和比特数

    在机器视觉系统中,反映每一个像元灰度质量的指标是动态范围,也是机器视觉系统要考虑的重要指标之一。动态范围和空间分辨率是机器视觉系统的两个主要指标。 灰度的动态范围在摄像头中的模拟视频部分用信噪比(Signal to Noise Ratio)SNR表示;而在摄像头或采集卡的A/

    2023年04月09日
    浏览(40)
  • Python 不同分辨率图像峰值信噪比[PSNR]

    PNNR:全称为“Peak Signal-to-Noise Ratio”,中文直译为峰值信噪比 前言 一、定义 二、Python代码 1.自定义 2.Tensorflow 总结 峰值信噪比是一种衡量图像质量的指标,描述的是最大值信号与背景噪音之间的关系。 一般来说,PSNR高于40dB说明图像质量极好(即非常接近原始图像);在

    2024年02月01日
    浏览(53)
  • 理解信噪比SNR,Eb/N0,Es/N0

    之前学习主要考虑的SNR和误码率,对Eb/N0和Es/N0不太了解,这次边记录边学习一下(希望随着自己的学习可以不断完善)。 信噪比,即信号功率与噪声功率的比值: S:信号功率,N:噪声功率。 Eb/N0:每个二进制bit能量与噪声功率谱密度的比值(比特信噪比)。 Es/N0:  每个符号

    2024年02月04日
    浏览(44)
  • 通信系统中基于matlab的BPSK信噪比检测算法及实现

    根据是否需要辅助数据,信噪比估计算法可以分为数据辅助类算法(Data aided, DA)和非数据辅助类算法(No Data aided, NDA)。DA估计算法准确性较高,但是需要提供先验信息,需要牺牲信道传输效率。NDA方法在传输数据信息的同时进行信噪比估计,不影响信息传输效率,适用范围较广

    2024年02月04日
    浏览(44)
  • python求不同分辨率图像的峰值信噪比,一文搞懂

    可以使用 Python 的 NumPy 和 OpenCV 库来实现这个任务。提前准备一张图片作为素材。 峰值信噪比(Peak Signal to Noise Ratio,PSNR)是衡量图像质量的常用指标,它表示图像中信号和噪声的比值。通常,较高的 PSNR 值表示图像质量较高。 PSNR 的公式如下: 其中, MAX 是图像的最大亮度

    2024年02月05日
    浏览(48)
  • 图像增强的两个评价指标:峰值信噪比PSNR和结构相似度SSIM

    图像增强的评价指标在像素层面上通常包含平均绝对误差(MAE)、均方误差法(MSE)、峰值信噪比(PSNR)以及结构相似度(SSIM)。目前在图像增强领域比较权威的客观评价标准为峰值信噪比(PSNR)和结构相似度(SSIM)。 注:这两个指标都需要由标准图做参考(不是原图)

    2024年02月05日
    浏览(46)
  • 适用于计算成像领域无参考图像的图像信噪比评价方法(SNR,PSNR,SSIM)(基础)

    适用于计算成像领域无参考图像的图像信噪比评价方法(基础) Image Signal-to-Noise-ratio evaluation method to reference-free images in the field of computitional imaging (basic). 注:英文可以不看 ,博主在练习英文而已,英文只是中文的翻译,可以直接看中文! 在许多计算成像领域中,我们没有办

    2024年02月03日
    浏览(44)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包