时间反转镜_matlab.zip

clc; clear all; %%%%%%%% 建立单频矩形脉冲信号 %%%%%%%% Time = 4.0000e-004; %脉冲宽度 Fsap = 4.e4; %采样频率 Num=256; %采样点数 t = (0:Num-1)*(1/Fsap); %时间参量 tri = zeros(1,Num); tra = t/Time; rect = zeros(1,Num); %矩形函数 for i = 1:Num tri(i) = tra(i)/max(tra); if (tri(i)>=-0.5)&&(tri(i)<=0.5) rect(i) = 1; end end f0 = 1.e4; %信号频率为10KHz Amp = 1; %信号幅值增益 signal = Amp*1/sqrt(Time)*rect.*exp(j*2*pi*f0*t); %单频矩形脉冲 figure(1); subplot(2,1,1); plot(t/1.e-3,real(signal));xlabel('时间/ms'); ylabel('S(t)/V');title('单频矩形脉冲基源信号'); NFFT = 2^nextpow2(Num); %FFT采样点数 f = Fsap/2*linspace(0,1,NFFT/2+1); Sw = fft(signal,NFFT)/Num; subplot(2,1,2); plot(f/1.e3,abs(Sw(1:NFFT/2+1)));xlabel('频率/kHz'); ylabel('S(w)/V');title('频域'); Num1 = 4; A1 = 1*rand(1,Num1); A2 = 0.5*rand(1,Num1); %A1，A2分别代表海洋各个信道的幅值衰减和时间延迟 yt1 = zeros(Num1,Num); yt2 = zeros(1,Num); for i = 1:Num1 yt1(i,:) = Amp*1/sqrt(Time)*rect.*exp(j*2*pi*f0*(t-A1(i)))*A2(i); yt2 = yt2+yt1(i,:); end noise = 10*normrnd(0,0.5,1,Num); %均值为0，方差为0.25（标准差为0.5）的高斯白噪声 yt = yt2+noise; % plot(t/1.e-3,real(yt));xlabel('时间/ms'); % ylabel('y(t)/V');title('单水听器接收信号'); YT = fft(yt,NFFT)/Num; % plot(f/1.e3,abs(YT(1:NFFT/2+1)));xlabel('频率/kHz'); % ylabel('Y(w)/V');title('频域'); yt3 = zeros(1,Num); for i = 1:Num1 yt1(i,:)=Amp*1/sqrt(Time)*rect.*exp(j*2*pi*f0*t)*(A2(i)^2); yt3 = yt3+yt1(i,:)+noise*A2(i); end figure(2); subplot(2,1,1); plot(t/1.e-3,real(yt3));xlabel('时间/ms'); ylabel('z(t)/V');title('单水听器时间反转处理'); NFFT = 2^nextpow2(Num); %FFT采样点数 f = Fsap/2*linspace(0,1,NFFT/2+1); YT3 = fft(yt3,NFFT)/Num; subplot(2,1,2); plot(f/1.e3,abs(YT3(1:NFFT/2+1)));xlabel('频率/kHz'); ylabel('Z(w)/V');title('频域'); [a1,b1] = xcorr(yt); % figure(3); % plot(b1,real(a1));title('第一次接收信号的自相关'); [a2,b2] = xcorr(yt3); figure(3); % subplot(2,1,1); plot(b2,real(a2));title('TRM处理后自相关'); Fs = Fsap; % subplot(2,1,2); % Hs = spectrum.welch; % psd(Hs,yt3,'Fs',Fs); NFFT1 = 2^nextpow2(length(b2)); f1 = Fsap/2*linspace(0,1,NFFT1/2+1); A2 = fft(a2,NFFT1)/(length(b2)); %自相关函数的傅立叶变换得到功率谱 % figure(4); % subplot(2,1,1); % plot(f1/1.e3,10*log10(abs(A2(1:NFFT1/2+1))));xlabel('频率/kHz'); % ylabel('P(w)/dB');title('功率谱'); window = hamming(30); noverlap = 20; [Pxx,f] = pwelch(yt3,window,noverlap,NFFT,Fs); % Pxx1 = 10*log10(Pxx); % subplot(2,1,2); % plot(f/1.e3,Pxx1);xlabel('频率/kHz'); % ylabel('P(w)/dB');title('功率谱')