-
Notifications
You must be signed in to change notification settings - Fork 3
/
figure_14.m
184 lines (145 loc) · 5.58 KB
/
figure_14.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
% Code to generate Figure 14 of the echo statistics tutorial
%
% This code plots the PDF and PFA of the echo magnitude due to various types
% of scatterers (one at a time) randomly distributed in the sensor beam:
% point scatterer, Rayleigh scatterer, randomly oriented smooth prolate
% spheroid, and randomly oriented rough prolate spheroid.
% 2D distribution of scatterers.
%
% Author: Kyungmin Baik | [email protected] | KRISS
clear all
close all
t0=clock;
N=1e4;
b=logspace(-5,0,N);
theta=[0:.001:90]*pi/180;
PDF=zeros(size(b));
%ka=2*pi;
ka=44.251078712;
% Beampattern PDF of a circular aperture
for m=1:length(b)-1
x=ka*sin(theta);
y=4*(besselj(1,x)).^2-b(m)*(x.^2);
thetaind=find(y(1:length(x)-1).*y(2:length(x))<0);
thetaval=(theta(thetaind)+theta(thetaind+1))/2;
PDF(m)=1/(2*pi*sqrt(b(m)))*sum(sin(thetaval)./(cos(thetaval).*abs(besselj(2,ka*sin(thetaval))))); % 2D
%PDF(m)=1/(4*sqrt(b(m)))*sum(((sin(thetaval)).^2)./(cos(thetaval).*abs(besselj(2,ka*sin(thetaval))))); % 3D
end
figure(1)
Nfac=sqrt(trapz(b,(b.^2).*PDF)); % normalization factor for fs
bran=b/Nfac;
PDFnorm=PDF/trapz(bran,PDF);
loglog(bran,PDFnorm,'k','LineWidth',1)
hold on
figure(2)
% Beampattern PFA of a circular aperture
PFA=zeros(size(PDFnorm)-1);
for m=1:N-1
PFA(m)=trapz(bran(m:N),PDFnorm(m:N));
end
loglog((bran(1:N-1)+bran(2:N))/2,PFA,'k','LineWidth',1)
hold on
ae=0.1; % radius of equal volume of sphere in m
fes=ae/2; % backscattering amplitude for equal volume of sphere.
% Combined with Rayleigh PDF
Nr=1e4;
Fray=logspace(-5,2,Nr);
Rran=Fray/fes;
PDFray=(Rran/fes).*exp(-0.5*(Rran.^2));
NT=1e4;
Fbeam=logspace(-5,0,NT);
%Ftotal=logspace(-20,log10(0.999999*Fmax),NT); % for smooth spheroid case
Ftotal=logspace(-9,2,NT);
% Combined with beampattern (Rayleigh scatterer)
PDFtotal=zeros(size(Ftotal));
for m=1:NT
Bran=fliplr(Ftotal(m)./Fray); % Range of beampattern pdf
%PDFtotal(m)=trapz(Bran,interp1(Ftotal,PDFbeam,Bran,'linear',0)./Bran.*fliplr(PDFray));
PDFtotal(m)=trapz(Bran,interp1(Fbeam,PDF,Bran,'linear',0)./Bran.*fliplr(PDFray));
end
Nfac=sqrt(trapz(Ftotal,(Ftotal.^2).*PDFtotal)); % normalization factor for fs
% PDF of Rayleigh scatterer combined with beampattern PDF
figure (1)
loglog(Ftotal/Nfac,PDFtotal/trapz(Ftotal/Nfac,PDFtotal),'b','LineWidth',2)
axis([1e-2 1e2 1e-6 1e4])
% PFA of Rayleigh scatterer combined with beampattern PDF
figure(2)
PFA=zeros(size(PDFtotal)-1);
for m=1:NT-1
PFA(m)=trapz(Ftotal(m:NT)/Nfac,PDFtotal(m:NT)/trapz(Ftotal/Nfac,PDFtotal));
end
loglog((Ftotal(1:NT-1)+Ftotal(2:NT))/(2*Nfac),PFA,'b','LineWidth',2)
axis([1e-2 1e2 1e-10 1e0])
PDFbeam=PDF;
Ep=0.1; % aspect ratio of spheroid
Fmax=fes*Ep^(-2/3);
Fmin=fes*Ep^(4/3);
Fran=logspace(log10(1.000001*Fmin),log10(0.999999*Fmax),N);
% PDF of spheroid in free-field
%PDF=((Ep^(2/3))*fes)./(2*(Fran.^(3/2))*sqrt(1-Ep^2).*sqrt(fes-(Ep^(2/3))*Fran)); % PDF for spheroid 3D
clear PDF
PDF=((Ep^(2/3))*fes)./(pi*Fran.*sqrt(fes-(Ep^(2/3))*Fran).*sqrt(Fran-(Ep^(4/3))*fes)); % 2D space
% This is for the calculation of smooth spheroid
PDFray=PDF;
Fray=Fran;
Fbeam=logspace(-5,0,NT);
Ftotal=logspace(-9,2,NT);
% Combined with beampattern (smooth spheroid)
PDFtotal=zeros(size(Ftotal));
clear Bran
for m=1:NT
Bran=fliplr(Ftotal(m)./Fray); % Range of beampattern pdf
%PDFtotal(m)=trapz(Bran,interp1(Ftotal,PDFbeam,Bran,'linear',0)./Bran.*fliplr(PDFray));
PDFtotal(m)=trapz(Bran,interp1(Fbeam,PDFbeam,Bran,'linear',0)./Bran.*fliplr(PDFray));
end
Nfac=sqrt(trapz(Ftotal,(Ftotal.^2).*PDFtotal)); % normalization factor for fs
% PDF of smooth spheroid combined with beampattern PDF
figure(1)
loglog(Ftotal/Nfac,PDFtotal/trapz(Ftotal/Nfac,PDFtotal),'g','LineWidth',2)
% PFA of smooth spheroid combined with beampattern PDF
figure(2)
clear PFA
PFA=zeros(size(PDFtotal)-1);
for m=1:NT-1
PFA(m)=trapz(Ftotal(m:NT)/Nfac,PDFtotal(m:NT)/trapz(Ftotal/Nfac,PDFtotal));
end
loglog((Ftotal(1:NT-1)+Ftotal(2:NT))/(2*Nfac),PFA,'g','LineWidth',2)
% Combined with Rayleigh PDF
clear Fray PDFray
Fray=logspace(-5,2,Nr);
PDFray=zeros(size(Fray));
for m=1:Nr
Rran=fliplr(Fray(m)./Fran); % Range of Rayleigh distribution
PDFray(m)=trapz(Rran,exp(-0.5*(Rran.^2)).*fliplr(PDF));
end
clear Fbeam Ftotal
Fbeam=logspace(-5,0,NT);
Ftotal=logspace(-9,2,NT);
% Combined with beampattern (Rough spheroid)
PDFtotal=zeros(size(Ftotal));
clear Bran
for m=1:NT
Bran=fliplr(Ftotal(m)./Fray); % Range of beampattern pdf
%PDFtotal(m)=trapz(Bran,interp1(Ftotal,PDFbeam,Bran,'linear',0)./Bran.*fliplr(PDFray));
PDFtotal(m)=trapz(Bran,interp1(Fbeam,PDFbeam,Bran,'linear',0)./Bran.*fliplr(PDFray));
end
Nfac=sqrt(trapz(Ftotal,(Ftotal.^2).*PDFtotal)); % normalization factor for fs
% PDF of Rough spheroid combined with beampattern PDF in 2D
figure(1)
loglog(Ftotal/Nfac,PDFtotal/trapz(Ftotal/Nfac,PDFtotal),'r','LineWidth',2)
% Add Rayleigh PDF
loglog(Fray,2*Fray.*exp(-1*(Fray.^2)),'k','LineWidth',2)
axis([1e-4 1e2 1e-6 1e5])
% PFA of Rough spheroid combined with beampattern PDF in 2D
figure(2)
clear PFA
PFA=zeros(size(PDFtotal)-1);
for m=1:NT-1
PFA(m)=trapz(Ftotal(m:NT)/Nfac,PDFtotal(m:NT)/trapz(Ftotal/Nfac,PDFtotal));
end
loglog((Ftotal(1:NT-1)+Ftotal(2:NT))/(2*Nfac),PFA,'r','LineWidth',2)
% Add Rayleigh PFA
loglog(Fray,exp(-1*(Fray.^2)),'k','LineWidth',2)
legend('Point scatterer','Rayleigh scatterer (1:1)','Smooth 10:1 prolate spheroid','Rough 10:1 prolate spheroid')
axis([1e-4 1e2 1e-3 1e1])
eval(['disp(''total ',num2str(etime(clock,t0)),' seconds elapsed'');'])