[READ-ONLY] Mirror of https://github.com/jmrplens/FDTDexamples. Modelos FDTD de aeroacústica simples y complejos para MATLAB 2019b o superior
acoustic
fdtd
matlab
Novedad
Jose Manuel Requena Plens
(Mar 16, 2020, 1:34 AM +0100)
424d5bfd
774f8723
+208
+208
FDTD_Columna_coord.m
+208
FDTD_Columna_coord.m
···
1
+
%% =========================================================================
2
+
%
3
+
% Modelo FDTD de una columna de altavoces direccionable.
4
+
%
5
+
% Se incluyen 3 puntos de recepci�n.
6
+
%
7
+
% Jose Manuel Requena Plens (16/03/2020)
8
+
%
9
+
% =========================================================================
10
+
clc, clear, close all;
11
+
% Version de la columna donde se maximiza la eficiencia en el punto de
12
+
% escucha optimizando inclinacion y curvatura
13
+
14
+
%% GENERAL
15
+
% Caracteristicas del medio
16
+
rho = 1.21; % Densidad del medio (kg/m3)
17
+
c = 341; % Velocidad de propagaci�n (m/s)
18
+
k = (c^2)*rho; % Modulo de compresibilidad
19
+
% Frecuencia de emision
20
+
freq = 4000; % (Hz)
21
+
% Tipo de excitaci�n ('pulse' o 'CW')
22
+
excit = 'CW';
23
+
% Dimensiones del mapa en metros
24
+
lx = 5; % m
25
+
ly = 3; % m
26
+
% Columna
27
+
N = 32; % Numero de elementos de la columna
28
+
d = 0.04; % Distancia entre elementos (de centro a centro) (m)
29
+
Origen = [1,2.5]; % Posicion de origen de la columna (m,m)
30
+
% Punto de escucha
31
+
Foco = [4,0.5];
32
+
% Posicion de los receptores
33
+
rec1 = [4,0.25]; % Posicion del receptor 1 (m,m)
34
+
rec2 = [1,0.25];% Posicion del receptor 2 (m,m)
35
+
rec3 = [lx-1,Origen(2)-(N-1)*d/2]; % Receptor a 0 grados
36
+
% Tama�o/Ancho de la PML en metros
37
+
lPML = 0.25; % m
38
+
% Par�metros de calculo
39
+
dh = .01; % Definici�n espacial (m)
40
+
dt = dh/341/2; % Definici�n temporal (s)
41
+
ts = 15; % Tiempo de simulaci�n (ms)
42
+
43
+
%% MATRICES
44
+
% Dimensiones de la matriz del mapa
45
+
nPML = round(lPML/dh);
46
+
nx = round(lx/dh)+nPML*2; % Por 2 ya que la PML esta arriba y abajo
47
+
ny = round(ly/dh)+nPML*2; % Por 2 ya que la PML esta a izquierda y derecha
48
+
% Crear matriz del mapa
49
+
p = zeros(nx,ny); % Matrices de presion
50
+
px = zeros(nx,ny);
51
+
py = zeros(nx,ny);
52
+
ux = zeros(nx+1,ny); % Matriz de velocidad de particula en x
53
+
uy = zeros(nx,ny+1); % Matriz de velocidad de particula en y
54
+
55
+
%% POSICIONES DE LOS ELEMENTOS
56
+
% Posiciones de los elementos en la matriz del mapa
57
+
posNy = round((Origen(2)-(0:N-1) * d)/dh) + nPML;
58
+
posNx = round(Origen(1)/dh) + nPML;
59
+
60
+
%% C�LCULO DE LOS RETARDOS TEMPORALES (dly)
61
+
% Calculo del angulo para emitir hacia el punto de escucha
62
+
theta = -rad2deg(atan((Origen(2)-d*N/2-Foco(2)+lPML)/(Foco(1)-Origen(1))));
63
+
% Inclinaci�n
64
+
x = sin(-theta*pi/180)*d; % Calculo de distancia simulada del segundo elemento
65
+
D = x * (0:N-1)'; % Vector de distancias para cada elemento (de arriba a abajo)
66
+
dlytilt = D/c; % Vector de tiempos/delays para cada elemento por la inclinacion
67
+
68
+
% Calculo de la distancia desde el origen de la columna (centro de la
69
+
% columna inclinada) hasta el punto de escucha para utilizarlo de radio
70
+
% para la focalizaci�n
71
+
rarc = -sqrt((Foco(1)-(Origen(1)-D(round(N/2))))^2+(Foco(2)-(Origen(2)-d*N/2))^2);
72
+
% Arco
73
+
thetaarc = linspace(-pi/2-deg2rad(-theta), pi/2-deg2rad(-theta), 1000); % Vector de angulos para una semicircunferencia
74
+
Ltot = (N-1)*d; % Longitud total de la columna
75
+
% Semicircunferencia
76
+
x = (rarc * cos(thetaarc) + Origen(1)-rarc);
77
+
y = (rarc * sin(thetaarc) + Origen(2)-Ltot/2);
78
+
% Distancias
79
+
Darc = zeros(N,1);
80
+
for m = 1:N
81
+
[~,idx]=min(abs(y-(Origen(2)-(m-1)*d)));
82
+
Darc(m) = Origen(1)-abs(x(idx));
83
+
end
84
+
dlyarc = Darc/c; % Vector de tiempos/delays producidos por arquear
85
+
86
+
% Vector de delays
87
+
dly = dlyarc+dlytilt;
88
+
89
+
% Representacion de las fuentes virtuales frente a la real
90
+
% posVirtual = zeros(N,2);
91
+
% for m = 1:N
92
+
% posVirtual(m,:) = [Origen(1)-D(m)-Darc(m),Origen(2)-(m-1)*d+lPML];
93
+
% end
94
+
% figb = figure('Color',[1,1,1]);
95
+
% plot(posVirtual(:,1),posVirtual(:,2),'ro', 'MarkerSize', 2,'MarkerFaceColor','r')
96
+
% hold on
97
+
% plot((ones(size(posNy))*posNx-nPML)*dh,posNy*dh,'ko', 'MarkerSize', 2,'MarkerFaceColor','k')
98
+
% axis equal,grid on
99
+
% axis([0,lx,0,ly])
100
+
% title('Posici�n de las fuentes')
101
+
% xlabel('X [m]'),ylabel('Y [m]')
102
+
% legend('Virtual','Real')
103
+
104
+
%% EXCITACI�N (la misma se�al con diferente retardo para cada elemento)
105
+
% El retardo se aplica en el vector 't'
106
+
lenT = ts*10^-3/dt; % Longitud del vector de tiempo
107
+
aa = freq/(sqrt(pi)/2)*4;
108
+
t=((1:lenT)/(1/dt)-4/aa); % Vector de tiempos
109
+
w = zeros(lenT,N);
110
+
for n = 1:N
111
+
% Funcion que define la excitacion, una para cada n, se incluye el
112
+
% retardo para cada elemento
113
+
switch excit
114
+
case 'pulse'
115
+
w(:,n) = -(exp(-aa^2*((t-dly(n)).^2)/2).*(aa^2*((t-dly(n)).^2)-1)); % Impulso / Ricker con retardo
116
+
case 'CW'
117
+
w(:,n) = cos(2*pi*freq*t)+1i*sin(2*pi*freq*t); % Seno complejo
118
+
% Retardo aplicado
119
+
w(:,n) = w(:,n).*exp(dly(n).*-1j*2*pi*freq);
120
+
end
121
+
end
122
+
123
+
%% PML
124
+
% Gradiente de impedancias desde el valor del medio hasta un porcentaje de
125
+
% este
126
+
gammamax = 0.5; % Maxima reducci�n de la impedancia / porcentaje
127
+
% PML Superior e Inferior
128
+
gammaux = zeros(nx+1,ny);
129
+
gammaux(1:nPML,:) = repmat(gammamax*((nPML:-1:1)'/nPML).^2,1,ny); % Izquierda
130
+
gammaux(1+end-nPML:end,:) = repmat(gammamax*((1:1:nPML)'/nPML).^2,1,ny); % Derecha
131
+
gammax = (gammaux(1:end-1,:)+gammaux(2:end,:))/2; % Conjunto
132
+
% PML Izquierda y Derecha
133
+
gammauy = zeros(nx,ny+1);
134
+
gammauy(:,1:nPML) = repmat(gammamax*((nPML:-1:1)/nPML).^2,nx,1); % Inferior
135
+
gammauy(:,1+end-nPML:end) = repmat(gammamax*((1:1:nPML)/nPML).^2,nx,1); % Superior
136
+
gammay = (gammauy(:,1:end-1)+gammauy(:,2:end))/2; % Conjunto
137
+
138
+
%% Inicializacion de los vectores de los puntos de escucha
139
+
h = zeros(1,lenT);
140
+
h2 = zeros(1,lenT);
141
+
h3 = zeros(1,lenT);
142
+
143
+
%% Condiciones de contorno
144
+
ux(1,:) = -p(1,:)/rho/c;
145
+
ux(end,:) = p(end,:)/rho/c;
146
+
uy(:,1) = -p(:,1)/rho/c;
147
+
uy(:,end) = p(:,end)/rho/c;
148
+
149
+
%% Calculo
150
+
fig1 = figure('Color',[1,1,1]);
151
+
152
+
colormap(jet(256))
153
+
for tt=1:lenT
154
+
% Presion
155
+
px = px.*(1-gammax)-k*dt/dh*diff(ux);
156
+
py = py.*(1-gammay)-k*dt/dh*diff(uy')';
157
+
158
+
% Excitaci�n
159
+
for n = 1:N
160
+
px(posNx,posNy(n)) = w(tt,n);
161
+
py(posNx,posNy(n)) = w(tt,n);
162
+
end
163
+
p = px+py;
164
+
165
+
% Velocidad
166
+
ux(2:nx,:) = ux(2:nx,:).*(1-gammaux(2:nx,:))-dt/rho/dh*diff(p);
167
+
uy(:,2:ny) = uy(:,2:ny).*(1-gammauy(:,2:ny))-dt/rho/dh*diff(p')';
168
+
169
+
% Respuesta al impulso
170
+
h(tt) = p(round(rec1(1)/dh)+nPML,round(rec1(2)/dh)+nPML);
171
+
h2(tt) = p(round(rec2(1)/dh)+nPML,round(rec2(2)/dh)+nPML);
172
+
h3(tt) = p(round(rec3(1)/dh)+nPML,round(rec3(2)/dh)+nPML);
173
+
174
+
% Representaci�n gr�fica
175
+
if tt/10==round(tt/10)
176
+
splmap = 10*log10(abs(p).^2/(2e-5)^2);
177
+
%splmap = splmap - max(splmap(:));
178
+
pcolor(((1:nx)-nPML)*dh,((1:ny)-nPML)*dh,splmap');
179
+
hold on;
180
+
plot([0,0,nx-2*nPML,nx-2*nPML,0]*dh,[0,ny-2*nPML,ny-2*nPML,0,0]*dh,'r--','linewidth',2)
181
+
plot(repmat(posNx-nPML,1,N)*dh,(posNy-nPML)*dh, 'wo', 'MarkerSize', 2,'MarkerFaceColor','w');
182
+
plot([rec1(1),rec2(1),rec3(1)],[rec1(2),rec2(2),rec3(2)], 'bo', 'MarkerSize', 4,'MarkerFaceColor','b');
183
+
text([rec1(1),rec2(1),rec3(1)]+0.05,[rec1(2),rec2(2),rec3(2)],{'Rec1','Rec2','Rec3'},'Color','white')
184
+
hold off
185
+
axis equal;axis([-nPML,nx-nPML,-nPML,ny-nPML]*dh)
186
+
set(gca,'Clim',[50 110]);shading flat,cc = colorbar;
187
+
title(['Tiempo = ' num2str(round((tt)*1000*dt)) ' ms']);
188
+
xlabel('X [m]'),ylabel('Y [m]'),cc.Label.String = 'SPL (dB)';
189
+
drawnow
190
+
end
191
+
end
192
+
% Respuesta al impulso / Punto de escucha
193
+
time = seconds(t);
194
+
maxx = max([abs(h),abs(h2),abs(h3)]);
195
+
h = h./maxx; h2 = h2./maxx; h3 = h3./maxx;
196
+
T = array2timetable([h',h2',h3'],...% Datos iniciales
197
+
'RowTimes',time); % Columna de tiempo
198
+
fig= figure('Color',[1,1,1]);
199
+
s = stackedplot(T,...
200
+
'Title','Se�al recibida',...
201
+
'DisplayLabels',{'Receptor 1', 'Receptor 2', 'Receptor 3'},...
202
+
'XLabel','Se�al recibida');
203
+
s.GridVisible = 1;
204
+
s.AxesProperties(1).YLimits = [-1.2,1.2];
205
+
s.AxesProperties(2).YLimits = [-1.2,1.2];
206
+
s.AxesProperties(3).YLimits = [-1.2,1.2];
207
+
fig.Children.FontName = 'Arial';
208
+
fig.Children.FontSize = 12;