%function Gausscos(choice) % can use to change the ways of imaging
% can be multiply choice:
% 1 :for Plane B-scan
% 2 : for Focused B-scan
% 3 : for Sector B-scan
close all
clear all
clc
choice = 1;
Fc=50000; % carrier frequency
Fs=1000000; % sample frequency
Ts=0.1; % sample time
Ns=Fs*Ts; % length of signal
t=0:1/Fs*10:0.003; % Generated signal domain
c_variance=0.01; % c in Gaussian signal
A0=2.5; %A0=3,full amplitude;A0<3 ;
a_coeff=1/((2*pi)^(1/2)*c_variance);
b_meanvalue=0.1; % mean value of Gaussian signal
c1=6300; % wave speed
% g1=a_coeff*exp(-(t-b_meanvalue).^2/(2*c_variance^2)); %Gauss signal
%
% g=g1; %original signal
% figure(1) % Generates a graphical window (1)
% subplot(2,1,1) % shape the figure into fluctuation two parts, will take on part
% plot(t,g) % In the above section plot original signal graphics
% gms=ammod(g,Fc,Fs,0,A0); % moduated signal
% subplot(2,1,2); % In the bottom section plot processed signals graphics
% plot(t, gms);
g1=exp(-(t-0.0015).^2/(2*0.00000001^2));
figure
plot(g1);
title('高斯信号');
s1=sin(2*pi*1000000*t);
figure
plot(s1); title('正弦信号');
g2=g1.*s1;
figure
plot(t, g2, 'k:');
title('调制信号');
gms = g2;
% x=ammod(gms,Fc,Fs,0,A0); % moduated signal
% figure;
% plot(x);
%
% figure(2);
% x1=fft(x,16192);
% mag=abs(x1);
% %phase
% f=(0:length(x1)-1)*Fs/length(x1);
% plot(f,mag)
% figure(3);
% According to the code above will GMS into a frequency domain to drawing
gms_pro=(fft(gms));
f=(0:length(gms_pro)-1)*Fs/length(gms_pro);
figure
plot(fftshift(gms_pro));
title('调制信号频谱');
% For every point and launch, calculating the distance receiver and reflection angles
pitch = 0.63; % element pitch
width = 0.53; % element width
space = pitch - width; % gap between element
num_of_ele = 4; % number of element, can be changed when use
ele_start = -num_of_ele/2*pitch+space/2+width/2; % coordinate of first element
htr_store = zeros(num_of_ele*num_of_ele,size(gms,2));
% save htr, num_of_ele*(num_of_ele-1) Represents the number of outer circular, Docked with the situation when receiver and transmitter in the same point
% size(gms,2) Each cycle represents the dimension of htr obtained
xt_store = zeros(num_of_ele*num_of_ele, 1); % Line index is alexandrine, used to store the abscissa denotes the launch
xr_store = xt_store; % Used to store the abscissa denotes the receiving end
Dtpr_store_L2 = xt_store; % Used to store launch end to reflect segment to receiving end distance
Dtpr_store_L3 = xt_store; %
flag = 1; % Indicators used to indicate the line index htr
lambda = 1.26; %wavelength
V = c1; %wave speed
for i=1:num_of_ele
for j=1:num_of_ele % the outer loop is the transmitter and receiving
Tn=ele_start + (i-1)*pitch; % transmitter horizontal ordinate, 0.00063 is step length, the same below
Rn=ele_start + (j-1)*pitch; % receiver
T=[Tn 0]; % transmitter coordinate
R=[Rn 0]; % receiver coordinate
xt_store(flag) = Tn; % save transmitter coordinate
xr_store(flag) = Rn; % save receiver coordinate
P=pitch*[0 -2-1e-2]; % reflected point in the medium
% P = c1*[0 -2-1e-2];
TR=R-T; %set TR,TP,PR as Vector
TP=P-T;
PR=R-P;
L1=sqrt(dot(TR,TR)); % Calculation TR length
L2=sqrt(dot(TP,TP)); % TP
L3=sqrt(dot(PR,PR)); % PR
Dtpr_store_L2(flag)=L2; %distance from T to P to R
Dtpr_store_L3(flag)=L3; %distance from T to P to R
at=acos(-P(2)/L2); % angles between TR,TP
ar=acos(-P(2)/L3);
a = width; % element width
ptx = sin(pi*a*sin(at)/lambda)/(pi*a*sin(at)/lambda); % function 4 from paper model ,transmitter
prx = sin(pi*a*sin(ar)/lambda)/(pi*a*sin(ar)/lambda); % receiver
Atr = A0/sqrt(L2*L3); %function 5
Htr = ptx*prx*Atr*gms_pro; %function 6
htr =(ifft(Htr));
htr = abs(hilbert(htr)); %function 7
% figure;plot(t+(L2+L3)/V, htr);
% htr = hilbert(htr);
% htr = abs(htr);
% figure(1);hold on;plot(t+(L2+L3)/V, htr,'k');
htr_store(flag, :) = htr; %save htr
flag = flag + 1; %update htr
% end
end
end
flag = 0; %reset flag
tmp_h = 0.01; %space
X = -1:tmp_h:1; % X,Z means medium space [-5,5]*[0,-10]
Z = 0-1e-2:-tmp_h:-2-1e-2; %move medium as -1e-2
[xx,zz] = meshgrid(X, Z); %X,Z meshgrid
Ixz = xx;Ixz(:) = 0; %I(x,z) reset
% figure
for i = 1:size(X, 2)
for j = 1:size(Z, 2) % loop to point in medium
P = [X(i) Z(j)]; % medium
D = 5; %function 9,10
for k = 1:num_of_ele*num_of_ele %loop of different transmitter and receiver
flag = flag + 1;%update the flag
T = [xt_store(flag), 0];%coordinate of transmitter
R = [xr_store(flag), 0];%coordinate of receiver
TP=P-T; %setTR,PR as Vector
PR=R-P;
L1=sqrt(dot(P,P)); %calculate the length between P and 0 point, used in case 3P L2=sqrt(dot(TP,TP)); %length of TP
L3=sqrt(dot(PR,PR)); %length of PR
s = 1.2; %function 10
[htr_max, htr_max_flag] = max(htr_store(flag,:));
switch choice
case 1 % case 1: Plane B-scan
% if (L2 < D/2) && ( L3 < D/2) %calculate domain
if ( L3 < D/2) %calculate domain
% dt = 2*Z(j)/c1 + (L2+L3-Dtpr_store(flag))/V; %time delay
% [t_real, t_real_flag] = min(abs(dt - t - Dtpr_store(flag)/V)); %most approximate time
dt = round((L2 + L3 - Dtpr_store_L2(flag) - Dtpr_store_L3(flag))/V*Fs/10);
t_real_flag = htr_max_flag + dt;
htr_real = htr_store(flag, t_real_flag);% htr value
Ixz(i, j) = Ixz(i, j) + htr_real;%update Ixz
end
case 2 % case 2: Focused B-scan
% if (L2 < D*s/2) && ( L3 < D*s/2)%calculate domain
if ( L3 < D*s/2) %calculate domain
% dt = (L2+L3)/c1 + (L2+L3-Dtpr_store(flag))/V;%time with delay
% [t_real, t_real_flag] = min(abs(dt - t - Dtpr_store(flag)/V));
%most approximate time
dt = round((L3 - Dtpr_store_L3(flag))/V*Fs/10);
t_real_flag = htr_max_flag + dt;
htr_real = htr_store(flag, t_real_flag);% htr value
Ixz(i, j) = Ixz(i, j) + htr_real; %update Ixz
end
case 3
% case 3: Sector B-scan
% dt = (2*L1 + (xt_store(flag) + xr_store(flag))*sin(atan2(X(i),Z(j))))/c1 + ...
% (L2+L3-Dtpr_store(flag))/V; %last result add on the time delay
%
评论5