A introdução de ruídos aditivos e multiplicativos nos modelos, tanto para o Pistão Circular Plano quanto para o Pistão Côncavo deve gerar resultados interessantes se verificados em relação a seus respectivos Erros Quadráticos médios. A escolha do tipo de geometria dos alvos pode determinar o tipo de modelo a ser adotado. Geometrias com células discretas ao invés de compactadas como nesse trabalho podem explorar suas resoluções axiais e radiais.
A utilização de mais transdutores transmitindo e recebendo sinais acústicos, tanto nas simulações quanto nos ensaios, pode fornecer mais informações acerca das propriedades dos alvos escolhidos.
Outros tipos de regularização podem ser estudados como por exemplo a regularização por variação total e as regularizações iterativas baseadas em gradiente conjugado.
REFERÊNCIAS BIBLIOGRÁFICAS
ARDITI, M.; FOSTER, F. S.; HUNT, J. W. Transient fields of concave annular arrays, Ultrasound Imaging, v. 3, 1981, p. 37-61.
BERNDT, H.; SCHNIEWIND, A. P.; JOHNSON, G. C. High-resolution ultrasonic imaging of wood, Wood Science and Technology, vol. 33, no. 3, June 1999, p. 185- 198.
BLAHUT, R. E. Theory of Remote Image Formation. Cambridge, MA: Cambridge University Press, 2004.
BOVIK, A. Handbook of Image and Video Processing. San Diego, CA: Academic Press, 2000.
BRIDAL, S. L. et al. Milestones on the road to higher resolution, quantitative, and functional ultrasonic imaging," in Proceedings of the IEEE, vol. 91, no. 10, October 2003, p. 1543-1561.
BRUSSEAU, E. et al. Fully automatic luminal contour segmentation in intracoronary ultra-sound imaging - a statistical approach, IEEE Transactions on Medical Imaging, vol. 23, no. 5, May 2004.
BURCKARD, C. B. Speckle in ultrasound B-mode scans, IEEE Transactions on Sonics and Ultrasonics, vol. 25, no. 1, January 1978, p. 1-6.
COLTON, D.; KRESS, R. Inverse Acoustic and Electromagnetic Scattering Theory, Springer, Berlin, 1992, p. 121, 289 (new edition: 1998, p. 133, 304).
ENGL, H. W.; HANKE, M.; NEUBAUER, A. Regularization of Inverse Problems. Boston, MA: Kluwer Academic Publishers, 1996.
GENE, H.; GOLUB, G. H.; MATT, U. v. Quadratically constrained least squares and quadratic problems, Numer. Math. 59. 1991.
GOLUB, G. H.; HEATH, M.; WAHBA, G. Generalized cross-validation as a method for choosing a good ridge parameter, Technometrics, vol. 21, 1979, p. 215-223. HANSEN, P. C. Analysis of discrete ill-posed problems by means of the L- curve,"SIAM Review, vol. 34, no. 4, December 1992, p. 561-580.
HANSEN, P. C. Regularization Tools: A Matlab package for analysis and solution of discrete ill-posed problems. Volume 6, Issue 1, 1994, p. 1-35.
HANSEN, P.C.; SEKII, T. The modified SVD method for regularization in general form, SIAM J. Sci. Statist. Comp. 13, 1992, p. 1142-1150.
HARRIS, G. R. Transient field of a baffled planar piston having an arbitrary vibration amplitude distribution. Journal of the Acoustical Society of America, v. 70, 1981, p. 186-204.
IVANOV, V. K. On linear problems which are not well-posed, Dokl. Akad. Nauk SSSR 145, 1962, p. 270–272 (Russian).
KETTERLING, J. A. Acoustic field of a wedge-shaped section of a spherical J. Acoust. Soc. Am., Vol. 114, No. 6, Pt. 1, Dec. 2003.
KINO, G. S. Acoustic Waves: Devices, Imaging, and Analog Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1987.
KINSLER, L. E. et al. Fundamentals of Acoustics. 4th ed. New York, NY: Wiley, 2000.
LAVARELLO, R. KAMALABADI, F.; O’BRIEN, WILLIAM D. A Regularized Inverse Approach to Ultrasonic Pulse-Echo Imaging. IEEE Transactions on Medical Imaging, v. 25, n. 6, 2006, p. 712-22.
LAWSON, C. L.; HANSON, R. J. Solving Least Squares Problems. Philadelphia, PA: Prentice-Hall, 1974.
LEE, H.; SCHUELER, C. F.; Wade, G. Fundamentals of digital ultrasonic imaging, IEEE Transactions on Sonics and Ultrasonics, vol. 31, no. 4, July 1984, p. 195-217. LOCKWOOD, J. C.; WILLETTE, J. G. High-speed method for computing the exact solution for the pressure variations in the nearfield of a baffled piston, J. Acoust. Sot. Amer. 53, 1973, p. 735-741.
MAHAUT, S. et al. Application of phased array techniques to coarse grain components inspection, Ultrasonics, vol. 42, nos. 1-9, April 2004, p. 791-796. MCLAREN, S.; WEIGHT, J.P. Transmit-receive mode responses from finite-sized targets in fluid media. Journal of the Acoustical Society of America, v. 82, n. 6, , 1987, p. 2102-2112.
MOROZOV, V. A. On the solution of functional equations by the method of regularization, Soviet Mathematics - Doklady, vol. 7, , 1966, p. 414-417.
MUELLER, J. L., SILTANEN, S. Linear and Nonlinear Inverse Problems Applications. SIAM, 2012.
OELZE, M. L. et al. Differentiation and characterization of rat mammary broadenomas and 4T1 mouse carcinomas using quantitative ultrasound imaging, IEEE Transactions on Medical Imaging, vol. 23, no. 6, , June 2004, p. 764- 771.
O'NEIL, H. T. Theory of focusing radiators, J. Acoust. Sot. Amer. 21, 516-526 .1949. PENTTINEN, P.; LUUKKALA, M. The impulse response and pressure nearfield of a curved ultrasonic radiator, J. Phys. D 9, 1976, p. 1547–1557.
PHILLIPS, D. L. A technique for the numerical solution of certain integral equations of the first kind, J. Assoc. Comput. Mach. 9, 1962, p. 84–97.
PIERCE, A. D. Acoustics: An Introduction to its Physical Principles and Applications. Woodbury, NY: Acoustical Society of America, 1989.
PRESS, W. H. et al. Numerical Recipes in C: The Art of Scientific Computing
(second edition). Cambridge University England EPress. 1992. p. 994. ISBN 0-521-
43108-5.
ROBINSON, D. E.; CHEN, C. F.; WILSON, L. S. Measurement of velocity of propagation from ultrasonic pulse echo data. Ultrasound Med. Biol., 9.1982, p. 413–420.
ROBINSON, D. E.; LEES, S.; BESS, L. Near field transient radiation patterns for circular pistons, IEEE Trans. Acoustics, Speech and Signal Processing ASSP-22, 1974, p. 395-403.
SCHICKERT, M. KRAUSE, M.; MULLER, W. Ultrasonic imaging of concrete elements using reconstruction by synthetic aperture focusing technique, Journal of Materials in Civil Engineering, vol. 15, no. 3, May/June 2003, p. 235-246.
SCHUELER, C. F.; LEE, H.; WADE, G. IEEE Trans. Sonics and Ultrasonics, SU-31, 1984, p.195-2 17.
STEPANISHEN, P. Transient radiation from pistons in an infinite planar baffle, J. Acoust. Soc. Am. 49, 1971, p. 1629–1683.
TIKHONOV, A. N. On the solution of ill-posed problems and the method of regularization, Dokl. Akad. Nauk SSSR 151. 1963, p. 501–504 (Russian).
VOGEL, C.R. Computational Methods for inverse problem, SIAM, 2002.
VOGT, M. et al. In vivo evaluation and imaging of skin elasticity applying high frequency (22 MHz) ultrasound, in IEEE Ultrasonics Symposium, , 2002, p. 1863- 1866.
WAGNER, R. F. et al. Statistics of speckle in ultrasound B-scans, IEEE Transactions on Sonics and Ultrasonics, vol. 30, no. 3, May 1983, p. 153-163.
WATKINS, D. S. Fundamentals of Matrix Computations. New York: Wiley, 2002. WEIGHT, J. P. Ultrasonic beam structures in fluid media. Journal of the Acoustical Society of America, v. 76, n. 4, 1984, p. 1184-1191.
WEIGHT, J. P.; HAYMAN, A. J. Observations of the propagation of very short ultrasonic pulses and their reflection by small targets. Journal of the Acoustical Society of America., v. 63, n. 2, 1978, p. 396-404.
ANEXOS
Anexo I - Rotinas de software do pistão plano teórico
1) Rotina montacampoacustico.m que monta o campo acústico:
%Monta a matriz S dos sistema S.r=g %
% Flavio Buiochi e Orlando Cirullo (jan/2013)
%%%%%% DADOS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
a=9.5e-3; %Raio do Pistão Plano Circular (m)
c=1500; %Velocidade de Propagação na água (m/s)
dens=1000; %Densidade (kg/m3)
freq= 1e6; %Frequencia central do transdutor (Hz)
dT=1/freq/5; %Amostragem (dT = 1/fa)
% Define a onda de excitaçao vn (experimental ou teórica)
periodo= 0:dT:(1/freq); %Período discretizado para um ciclo senoidal
[O]=Onda_a(5,9,2.4,1/freq,dT) vn = O.onda; %sinal de excitação
% Pontos a serem varridos
x = [-16*1e-3:0.25e-3:16*1e-3]; z = [50*1e-3:0.25e-3:70*1e-3];
% Posições de deslocamento do transdutor (emissor/receptor)
u = linspace(-10e-3,10e-3,21);
ti=73e-6; %66.675e-6; %tempo inicial(us)
tf=95e-6; %115.000e-6; %109.625e-6; %tempo final (us)
tic
[X,Z] = meshgrid(x,z); %define o conjunto de pontos na malha retangular
% empilhando todos os pontos do ROI
x=reshape(X,[prod(size(X)) 1]); %transforma em vetor as coords. x
z=reshape(Z,[prod(size(Z)) 1]); %transforma em vetor as coords. z
[Nz, Nx]=size(X); %numero de pontos nas direções x e z do ROI
N = length(x); %numero de pontos da malha
for i=1:N
i, x(i), z(i)
[phi, t0] = Potencial(a,c,x(i),z(i),dT);
% Cálculo da Pressão Impulsiva Normalizada
h=dens*diff([0 phi 0])/dens/c; pe = conv(vn,h);
if 0
te = t0*c*1e3 + (0:length(pe)-1)'*dT*c*1e3; %tempo (c.t (mm))
figure(1)
plot(te,pe,'r'); xlabel('c.t (mm)');
ylabel('Amplitude normalizada em volts'); title('Pulso transmitido ');
grid pause end pulsos(i).pe = pe; pulsos(i).t0 = t0;
% Cálculo da Amplitude de Pressão pico-a-pico p1=min(pe); p2=max(pe); campo(i)=abs(p1)+abs(p2); end if 0 for i=1:N i, x(i), z(i) pe = pulsos(i).pe;
te = t0*c*1e3 + (0:length(pe)-1)'*dT*c*1e3; %tempo (c.t (mm))
figure(1)
plot(te,pe,'r'); xlabel('c.t (mm)');
ylabel('Amplitude normalizada em volts'); title('Pulso transmitido ');
grid pause end
end
P=reshape(campo,size(X)); %retorna os dados do vetor na forma matricial
%Plotagem do campo acústico
figure (1)
surf(X,Z,P) , shading interp
xlabel('x [mm]'); ylabel('z[mm]');
title('Campo Acústico - Imagem Original'); colorbar
shading flat
set(gcf, 'renderer', 'zbuffer') figure (2)
surf(X,Z,P) , shading interp
xlabel('x [mm]'); ylabel('z[mm]');
title('Campo Acústico - Imagem Original'); colorbar
shading flat
set(gcf, 'renderer', 'zbuffer') view(2)
axis equal
toc
2) Rotina monta_g_teorico.m que monta vetor g (simula os dados experimentais)
%
% Flavio Buiochi e Orlando Cirullo (Fev/2014)
Narq = 21; %número de sinais A-scan (experimental)
% Parâmetros de Tempo
ti=73e-6; %tempo inicial(us)
tf=95e-6; %tempo final (us)
dt=0.12500e-6; %discretização no tempo
tp=ti:dt:tf;
g = [];
for col = 1:Narq %varrendo as posições [u1, u2... uNarq] do
transdutor
load (['sinal_' num2str(col)]) % carrega sinais A-scan
method = 'linear';
yi = interp1(te,Pe,tp,method); %gera o vetor g de t1=ti a tp<=tf com dt
g = [g; yi']; % Empilha em vetor coluna
end
3) Rotina monta_S_teorico.m que Monta a matriz S dos sistema S.r=g
%
% Flavio Buiochi e Orlando Cirullo (jan/2013)
%%%%%% DADOS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=10e-3; %Raio do Pistão Plano Circular (m)
cc =1499.25%Velocidade de Propagação na água (m/s)
dens=1000; %Densidade (kg/m3)
freq=1e6; %Frequencia central do transdutor (Hz)
dT=1/freq/20; %Amostragem (dT = 1/fa)
% Define a onda de excitaçao vn (experimental ou teórica)
periodo= 0:dT:(1/freq); %Período discretizado - um ciclo senoidal
[O]=Onda_a(5,9,2.4,1/freq,dT) vn = O.onda; %sinal de excitação
% % Pontos do ROI a serem varridos
x = [-16*1e-3:1e-3:16*1e-3]; z = [50*1e-3:1e-3:70*1e-3];
% % Posições de deslocamento do transdutor (emissor/receptor)
u=linspace(-10e-3,10e-3,21); %u=[-10e-3, 0, 10e-3];
%% Parâmetros de Tempo
ti=73e-6; %tempo inicial(us)
tf=95e-6; %tempo final (us)
dt=0.12500e-6; %discretização no tempo
tp=ti:dt:tf;
tic
[X,Z] = meshgrid(x,z); %define o conjunto de pontos na malha retangular
% empilhando todos os pontos do ROI
x=reshape(X,[prod(size(X)) 1]); %transforma em vetor as coords. x
z=reshape(Z,[prod(size(Z)) 1]); %transforma em vetor as coords. z
[Nz, Nx]=size(X); %numero de pontos nas direções x e z do ROI
N = length(x); %numero de pontos da malha
NTx = length(u); %numero de posições do transdutor
PE = [], S=[];
for q=1:NTx for i=1:N
[i, x(i), z(i)]
% Cálculo do potencial impulsivo
[phi, t0] = Potencial(a,cc,x(i)-u(q),z(i),dT); t0 = 2*t0;
%figure(1), plot(phi)
% Cálculo da Pressão Impulsiva Normalizada
h=dens*diff([0 phi 0])/dens/cc; % figure(2), plot(h) % Cálculo do pulso-eco pe = conv(vn,(conv(h,h))); t = t0 + dT*(0:1:length(pe)-1); if 0 t = t0 + dT*(0:1:length(pe)-1); figure(3)
plot(t*1e3*c,pe,'.r'); xlabel('c.t (mm)');
ylabel('Amplitude normalizada'); title('Pulso-Eco (vn*h*h)'); grid % pause end method = 'linear'; pei = interp1(t,pe,tp,method); peiaux = isnan(pei); I = find(peiaux == 1); pei(I) = 0; if 0 figure(4) hold on plot(tp*1e3*c,pei,'.k'); xlabel('c.t (mm)');
ylabel('Amplitude normalizada');
title('Pulso-Eco INTERPOLADO (vn*h*h)'); hold off pause end %limitar os vetores de ti a tf PE = [PE pei']; if 0 for ii=1:N; figure(5) hold on plot(PE(:,ii));
xlabel('c.t (mm)'); ylabel('Amplitude PE'); title('Pulso-Eco'); hold off pause; end end end S = [S; PE]; PE = []; end save matriz_S S toc
4) Rotina PE6.m que calcula o pulso-eco de um conjunto de refletores pontuais com coeficiente de reflexão R(x,z) para diversas posições do transdutor (emissor/receptor)
%
% Flavio Buiochi e Orlando Cirullo (julho/2014)
geom = 1;
d=20e-3; %Diâmetro do Pistão Circular (m)
a=d/2; %Raio do Pistão Plano Circular
c=1500; %Velocidade de Propagação na água (m/s)
dens=1000; %Densidade (kg/m3)
freq=1e6; %Frequencia central do transdutor (Hz)
T=1/freq;
dT=T/40; %Amostragem (dT = 1/fa)
periodo= 0:dT:(1/freq); %Período discretizado- um ciclo senoidal
[O]=Onda_a(5,9,2.4,1/freq,dT) vn = O.onda; %sinal de excitação
%% Pontos a serem varridos (definindo ROI) % Valores de u, x, z x = [-16*1e-3:0.25e-3:16*1e-3]; z = [50*1e-3:0.25e-3:70*1e-3]; u = linspace(-10e-3,10e-3,21); tic
[X,Z] = meshgrid(x,z); %define o conjunto de pontos numa malha retangular
R = zeros(size(X)); %Cria matriz de zeros para R de coeficientes de reflexão
%% GEOMETRIA PONTOS
% R(5,7)=1; R(1,7)=1; R(1,4)=0; R(5,4)=0;
%% GEOMETRIA RETA (geometria 0)
if geom == 0
[m n]=find(X>=-0.005 & X<=0.005 & Z>0.0600 & Z<0.0602); %[m n]=find(X>=-0.005 & X<=0.005 & Z==0.0600);
R(m,n)=1; end
%% GEOMETRIA 1: 2 segmentos de 5mm espaçados de 12mm
if geom ==1
[m n]=find(X>=-0.011 & X<=-0.006 & Z>0.0600 & Z<0.0602); R(m,n)=1;
[m n]=find(X>=0.006 & X<=0.011 & Z>0.0600 & Z<0.0602); R(m,n)=1;
end
%% GEOMETRIA 2: 2 segmentos de 5mm em profundidades diferentes
[m n]=find(X>=-0.011 & X<=-0.006 & Z==0.055); R(m,n)=1;
[m n]=find(X>=0.006 & X<=0.011 & Z==0.065); R(m,n)=1; % R(m,n)=0.5; end %% Plota geometria figure(1) mesh(X,Z,R), xlabel('x [mm]'); ylabel('z[mm]');
title('Imagem Original');
%axis equal, %view(2);
x=reshape(X,[prod(size(X)) 1]); %transforma em vetor as coords. x
z=reshape(Z,[prod(size(Z)) 1]); %transforma em vetor as coords. z
r=reshape(R,[prod(size(R)) 1]); %transforma em vetor os coefs de reflexão
N = length(x); %numero de pontos da malha (ROI)
NTx = length(u); %numero de posições do transdutor
Lmax = 0; %inicialização - comprimento máximo do vetor Pe
for q=1:NTx
u(q) %mostra as posições do transdutor
Vdmax = sqrt((a + abs(x-u(q))).^2 + z.^2); %calcula o vetor das distancias
Idmax = find(Vdmax == max(Vdmax)); % Define os Indices de distancia Vdmax iguais ao seu max valor
zmax = Vdmax(Idmax(1)); %Define o maximo valor de z ...
Itmin = round(2*min(z)/(c*dT)); % Define-se os Indices de tempo mínimo de ...
Itmax = round(2*zmax/(c*dT)); % Define-se os Indices de tempo máximo de ...
t0min = Itmin*dT; %instante inicial que a matriz sparse pode ter sinal.
NdT = Itmax - Itmin + length(vn) + 2;% Diferença de índices de tempo ou número de pontos do vetor tempo.
Ps = sparse(NdT,1); % Cria matriz sparse de dimensões (m,n)=(NdT,1)
for i=1:N
if r(i) ~= 0, %evita calcular o potencial para os pontos com R=0.
%i, x(i), z(i)
% Cálculo do potencial impulsivo
[phi, t0] = Potencial(a,c,x(i)-u(q),z(i),dT); %figure(1), plot(phi)
% Cálculo da Pressão Impulsiva Normalizada
h=dens*diff([0 phi 0])/dens/c; % figure(2), plot(h)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if 0, %q==2 t = t0 + dT*(0:1:length(h)-1); figure(2) subplot(211); plot(t*1e3*c,h,'r'); xlabel('c.t (mm)');
ylabel('Pressão Impulsiva Normalizada por
dens*c');
title('Pressão Impulsiva (p_i)'); grid E = 1*conv(vn,(conv(h,h))); t = 2*t0 + dT*(0:1:length(E)-1); subplot(212); plot(t*1.e3*c,E,'r'); xlabel('c.t (mm)');
ylabel('Amplitude normalizada'); title('Pulso-Eco (vn*p_i*p_i)'); grid
pause end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Pulso-eco - com sparse
hh = r(i)*conv(h,h); It0 = ceil(2*t0/dT); Ps = Ps + sparse((It0-Itmin+1):(It0- Itmin+length(hh)),1,hh,NdT,1); end end Pe = conv(vn,full(Ps));
te = t0min + (0:length(Pe)-1)'*dT; %tempo para o vetor sparse (s)
figure(3)
plot(te*c*1e3,Pe,'r'); xlabel('c.t (mm)');
ylabel('Amplitude normalizada em volts'); title('Pulso-Eco ');
grid
savefile = ['sinal_' num2str(q) '.mat'];
save(savefile,'te','Pe') %Salva sinais A-scan (Pulso-eco)
if length(Pe)>Lmax, %buscando o comprimento máximo
Lmax = length(Pe); end
end
% coloca zeros no final dos vetores Pe para ter mesmo % comprimento e atualiza os vetores te
for q = 1:NTx,
Peaux = zeros(1,Lmax);
load (['sinal_' num2str(q)]),
[te(1) te(length(te))]*1e6, %mostra o tempo inicial e final em us
Peaux(1:length(Pe)) = Pe; Pe = Peaux;
% Pe = awgn(Pe,10); %adicionando ruido branco gaussiano de 10dB no sinal Pe
te = t0min + (0:length(Pe)-1)'*dT;
savefile = ['sinal_' num2str(q) '.mat'];
save(savefile, 'te', 'Pe') % Salva sinais A-scan (Pulso-eco) c/ novos comprimentos
end
5) Rotina Potencial.m que calcula o potencial impulsivo do pistão plano
%
% Flavio Buiochi e Orlando Cirullo (jan/2012)
function [phi,t0] = Potencial(a,c,x,z,dT) x = abs(x);
z = abs(z);
%Calculo dos tempos limites
t0=z/c;
t1=sqrt((a-x)^2+(z^2))/c; t2=sqrt((a+x)^2+(z^2))/c;
t = t0:dT:t2; %Calculo do Vetor Tempo
Omega=zeros(1,length(t)); %inicializa, com zeros, o vetor Ômega(ct)
% Cálculo dos Arcos de Circunferência Omega(ct)
if x<a % dentro da superfície do pistão
%disp('entrei em x<a')
It1 = find(t0 <= t & t <= t1); It2 = find(t1 < t & t <= t2);
Omega(It1)=2*pi*ones(1,length(It1));% se t0<=t<=t1
Omega(It2)=2*acos((c^2*t(It2).^2-z^2+x^2- a^2)./(2*x*sqrt(c^2*t(It2).^2-z^2)));% se t1<t<=t2
elseif x==a % na borda do pistão
%disp('entrei em x=a')
It3 = find(t0==t & t == t1); It4 = find(t1 < t & t <= t2);
Omega(It3)=pi*ones(1,length(It3)); % se t=t0=t1
Omega(It4)=2*acos(sqrt(c^2*t(It4).^2-z^2)/(2*a));% se t1<t<=t2
elseif x>a % fora da superfície do pistão
%disp('entrei em x>a')
It5 = find(t1 <= t & t <= t2);
Omega(It5)=2*acos((c^2*t(It5).^2-z^2+x^2- a^2)./(2*x*sqrt(c^2*t(It5).^2-z^2))); % se t1<=t<=t2
end
phi=c*Omega./(2*pi); % Potencial de Velocidade
6) Rotina variavelocidadepropagacao.m que calcula a Regularização de Tikhonov para diversas velocidades de propagação
cc = 1500;
c_correto = 1;
processa; %usado para calcular o paramNorm quando c_correto =1
if c_correto %calcula o valor máximo para a solução Rideal usaado
no progr. 'Metricas.m' paramNorm = max(max(abs(RR))); c_correto=0; end cc = 1499 % c =1500.05; %1499.999; %delta_cc=0.00005; delta_cc=0.25; NN = 41;
clear Cplot, clear MSEplot
for j = 1:NN cc processa; Metricas; %pause Cplot(j)=cc; MSEplot(j)=MSE; cc = cc + delta_cc; end figure (13)
plot (Cplot,MSEplot,'-o') %(Cplot(1:NN),MSEplot(1:NN)) %como armazenar o MSE para cada c
xlabel('c [m/s]'); ylabel('MSE[%]');
7) Rotina SolucaoTikhonov.m que calcula a Solução do sistema S*r =g com g teórico (simulado)
% R é a matriz solução do vetor empilhado r % S é a matriz pulso-eco
% g é o vetor dos dados experimentais %
% Deve-se, primeiro carregar:
%load('C:\Users\Cirullo\Documents\MATLAB\Buiochi\Últimos
arquivos\monta_Sg_Tarugo0a19\matriz_S.mat') %??? carrega monta_S.m;
%
load(['C:\Users\Cirullo\Documents\MATLAB\Cirullo\Reg_monta_Sg_ Tarugo0a19\vetor_g_Teorico']) %??? carrega monta_g_Teorico.m; % ou
% load (['C:\Users\Cirullo\Documents\MATLAB\Buiochi\Últimos arquivos\monta_Sg_Tarugo0a19\sinal_' ])%??? carrega sinais A- scan teóricos
% solução não-regularizada por mínimos quadrados desde que(S'S)^?1 exista.
% Least Square % r = S\g;
% r =inv(S'*S)*S'*g;
% solução regularizada por mínimos quadrados desde que(S'S)^-1 exista.
% Gradiente Matrix L
for i = 1: length(S'*S)-1 L(i, i+1) = -1;
end
alfa = 1;
r_Tikh = inv(S'*S + alfa^2*L'*L)*S'*g;
% r_Tikh = inv(S'*S + alfa^2*L'*L)*S'*g/dens/c;
% r_Tikh = (pinv(S)*g); % Utlizando pseudoinversa de S
%% Método de Regularização de Tikhonov
r_Tikh = inv(S'*S + alfa^2*L'*L)*S'*g;
% r_Tikh = inv(S'*S + alfa^2*L'*L)*S'*g/dens/c;
% r_Tikh = (pinv(S)*g); % Utlizando pseudoinversa de S
%%Métodos de Regularização usando o Gradiente Conjugado % r_cg = conjgrad(S'*S,g); % r_cg = conjgrad(A,b,tol); % r_pcg = pcg(S'*S,g); % r_pcg = pcg(A,b,tol)
% r_bicg(S'*S,g); % r_bicg(A,b,tol); %%
%retorna os dados do vetor na forma matricial
XX = reshape(x,size(X)); ZZ = reshape(z,size(Z));
%R = reshape(r,size(X));
RR = reshape(r_Tikh,size(X));
% Plotando a Somatória dos dados na forma matricial(sem Normalização)
figure(7)
surf(XX,ZZ,RR) , shading interp
colorbar shading flat
xlabel('x[mm]'); ylabel('z[mm]');
title('Tikhonov - Transdutor Circular Plano ');
%set(gcf, 'renderer', 'zbuffer') %view(2)