main_zfft.m

利用频率细化方位估计,水声信号处理中的一项技术

%矢量水听器频率细化方位估计——利用zfft方法

clc
clear all
close all
%参数初始化
parameter_initial;

for T=1:time_length
    
    %矢量水听器信号模拟(远场条件，球面波衰减)
    [p,vx,vy]=vector_locomotion(x_a,y_a,target_x,target_y,target_vx,target_vy,T,r_0,len,fs,f_0,An,SNR);
    %改变fft点数
    nfft=fs;
    win=hanning(nfft); %窗函数
    ind_up=round(frequence_up/(fs/nfft));	%频率下限下标
    ind_down=round(frequence_down/(fs/nfft));	%频率上限下标
    %计算功率谱
	[sp,f]=psd(p,nfft,fs,win);
	[svx,f]=psd(vx,nfft,fs,win);
	[svy,f]=psd(vy,nfft,fs,win);
	[spvx,f]=csd(p,vx,nfft,fs,win);
	[spvy,f]=csd(p,vy,nfft,fs,win);
    %声强方法方位估计
    [amplitude_theta,theta,result]=azi_estimate_max(sp,svx,svy,spvx,spvy,nfft,fs,ind_up,ind_down);
    theta1_temp(T)=result;
    
    %改变fft点数
    nfft=fs*10;
    win=hanning(nfft); %窗函数
    ind_up=round(frequence_up/(fs/nfft));	%频率下限下标
    ind_down=round(frequence_down/(fs/nfft));	%频率上限下标
    %计算功率谱
	[sp,f]=psd(p,nfft,fs,win);
	[svx,f]=psd(vx,nfft,fs,win);
	[svy,f]=psd(vy,nfft,fs,win);
	[spvx,f]=csd(p,vx,nfft,fs,win);
	[spvy,f]=csd(p,vy,nfft,fs,win);
    %声强方法方位估计
    [amplitude_theta,theta,result]=azi_estimate_max(sp,svx,svy,spvx,spvy,nfft,fs,ind_up,ind_down);
    theta2_temp(T)=result;

    %频率细化
    nfft=fs*10;
	f_step=0.1;
	f_down=98;
	f_up=102;
	[fft_p,zf]=f_zfft(p,nfft,fs,f_step,f_down,f_up);
	[fft_vx,zf]=f_zfft(vx,nfft,fs,f_step,f_down,f_up);
	[fft_vy,zf]=f_zfft(vy,nfft,fs,f_step,f_down,f_up);
	
	%计算功率谱
	%自功率谱
	sp_zfft=abs(fft_p).^2/2;
	svx_zfft=abs(fft_vx).^2/2;
	svy_zfft=abs(fft_vy).^2/2;
	%互功率谱
	spvx_zfft= fft_p.*conj(fft_vx)/2;
	spvy_zfft= fft_p.*conj(fft_vy)/2;
    theta=atan2(real(spvy_zfft),real(spvx_zfft))*180/pi;				%求出各个频率指定的方位
    amplitude_theta=real(sp_zfft);
    ind_max=find(amplitude_theta==max(amplitude_theta));
    theta3_temp(T)=theta(ind_max);      %能量最大的方位信息


T

end;


error_1=abs(theta1_temp-45);
error_2=abs(theta2_temp-45);
error_3=abs(theta3_temp-45);
% error_4=abs(theta4_temp-45);
mean(error_1)
mean(error_2)
mean(error_3)
% mean(error_4)
% 
% std(theta1_temp-45);
% std(theta2_temp-45);
% std(theta3_temp-45);
% std(theta4_temp-45);
% 
% 
figure(1)
plot((1:time_length),theta1_temp,'b',...
    (1:time_length),theta2_temp,'r');
title('角度检测结果');
xlabel('x轴方位（m）');
ylabel('y轴方位（m）');

figure(2)
plot((1:time_length),error_1,'b',...
    (1:time_length),error_2,'r');
title('角度检测误差');
xlabel('x轴方位（m）');
ylabel('y轴方位（m）');