removing for loop by using 3d matrix

Hi. I wrote a code. I want to calculate pressure(p) of waves which produces by 1d array of elements. p is 2d variable. In the code I used 2d matrix and 2 loop. now I want to make a 3d matrix to remove one of the loops. but the problem that I have is when I creat 3d matrix p also gets a 3d matrix while it must be a 2d matrix.

clc;
close all;
clear all;
a=1;
c=1.5; % mm/us
T=5; dt=0.8; % us
fc=1;BW=3; % MHz
R=10;
xf=20;zf=20;
td=0;
dx=0.3; dz=0.3;pitch=0.3;
% x0=10; z0=0;
f=@(x) exp(-x.^2*BW^2);
g=@(x) sin(2*pi*fc*x);
H=@(x) f(x).*g(x);

Nz=120;
Nx=120;
x=[1:Nx]*dx;
z=[1:Nz]*dz;
xl=[1:Nx]*dx;
[z3,x3,xxl] = ndgrid(z,x,xl); 

for t=0:dt:T
            d3=sqrt((x3-xxl).^2+(z3).^2); 
            df=sqrt((xxl-xf).^2+(zf)^2);
            td=(R-df)./c;
            tt_hat=(d3./c)+td;
            p = 1./sqrt(d3).*H(t-tt_hat); %        
    imagesc(p(:,:,100)); 
%     colormap gray; colorbar;
%     set(gca,'clim');
%     title(['Time= ',num2str(t)]);
    pause(0.05);   
end
 

[Jan]

"when I creat 3d matrix p also gets a 3d matrix while it must be a 2d matrix"

This does not allow to understand, what you want to achieve. Seeing the code and knowing, that you want something else, is not enough to suggest a solution.

Please explain with details, what you want to change.

[Parisa Salemi]

Sorry. Thank you Jan. I want to produce foucsed beam. For this I had 4 loops at first. But by using 2Dmatrix I was able to reduce to two. Now I want to use 4d matrix to write my code without using for loops. I will also sent my previous code. The code I have to write should have the same results as the code below. I hope this time my explaination be clear.

clc;
close all, clear all;
a=1;
c=1.5; % mm/us
T=15; dt=0.8; % us
fc=1;BW=3; % MHz
R=10;
xf=20;zf=20;
td=0;
dx=0.3; dz=0.3;pitch=0.3;
x0=10; z0=0;
f=@(x) exp(-x.^2*BW^2);
g=@(x) sin(2*pi*fc*x);
H=@(x) f(x).*g(x);

Nz=120;
Nx=120;
x=[1:Nx]*dx;
z=[1:Nz]*dz;
%u=zeros(Nz,Nx);
[zz,xx] = ndgrid(z,x); 

for t= 0:dt:T
    sum1=zeros(Nz,Nx);
    for ind= 1:120
            x0=ind*pitch;
            dd=sqrt((xx-x0).^2+(zz-z0).^2); 
            df=sqrt((x0-xf)^2+(z0-zf)^2);
            td=(R-df)/c;
            tt_hat=dd/c+td;
            p = 1./sqrt(dd).*H(t-tt_hat); % 
            %p = 1./sqrt(dd).*(exp(-(t-tt_hat).^2*BW^2).*sin(2*pi*fc*(t-tt_hat))); 
            sum1=sum1+p;
    end
    cla;imagesc(sum1.^2); 
%     colormap gray; colorbar;
%     set(gca,'clim');
%     title(['Time= ',num2str(t)]);
    pause(0.05);
end

[Parisa Salemi]

I could to remove one of the loops. But now the problem that I have is, I do not know how can I remove the time loop in this code.

clc;
close all;
clear all;
a=1;
c=1.5; % mm/us
T=20; dt=0.8; % us
fc=1;BW=3; % MHz
R=10;
xf=20;zf=20;
td=0;
dx=0.3; dz=0.3;pitch=0.3;
% x0=10; z0=0;
f=@(x) exp(-x.^2*BW^2);
g=@(x) sin(2*pi*fc*x);
H=@(x) f(x).*g(x);
Nz=120;
Nx=120;
x=[1:Nx]*dx;
z=[1:Nz]*dz;
xl=[1:Nx]*dx;
[z3,x3,xxl] = ndgrid(z,x,xl);
xf1 = ones(size(xxl))*xf;
zf1 = ones(size(xxl))*zf;
    d3=sqrt((x3-xxl).^2+(z3).^2);
    df=sqrt((xxl-xf1).^2+(zf1).^2);
    td=(R-df)./c;
    tt_hat=(d3./c)+td;
for t=0:dt:T
    p = 1./sqrt(d3).*H(t-tt_hat);
    p1 = sum(p,3).^2;
    imagesc(p1);
    pause(0.05);
end

Answer 1

[Jan]

What is your purpose of vectorizing the code? The processing time is dominated by imagesc and pause here. I assume the loop are faster than the vectorizde code, which produces large intermediate arrays. Your original loop takes 0.86 seconds, the code in the comment above 1.33 seconds, if imagesc and pause are removed.

You can add an additional dimension:

dim = [1, size(tt_hat)];
p   = H((0:dt:T).' - reshape(tt_hat, dim)) ./ reshape(sqrt(d3), dim);
p1  = sum(p, 4).^2;

But I do not see the advantage compared to clean loops.

If efficiency matters, remember that anonymous functions are expensive. So reduce them to the minimum:

% This:
f=@(x) exp(-x.^2*BW^2);
g=@(x) sin(2*pi*fc*x);
H=@(x) f(x).*g(x);

% slows down the total processing time by 10% compare to this:
H = @(x) exp(-x.^2 * BW^2) .* sin(2 * pi * fc * x);

Note: Please format your code in the forum to improve the readability. I've done this for you today.

The brute clearing header "clc;close all;clear all;" is extremly inefficient: the clear('all') removes all functions from the memory and forces Matlab to reload them from the slow disk. This wastes a lot of time and offers no advantage. Prefer using functions instead to keep your workspace clean.

[Parisa Salemi]

So many thanks Jan for your great expolination. Sorry I can not understand what you mean about vectorizing the code. But my goal is creating an elament of array which produce wave and then see how the wave propagate in the x,y area. As I said before I did this by using for loops in my code. Now I am looking for any other way to do iteration in code without using for loop. My supervisour suggested to use 4Dmatrix.

Thanks for the code but in my case I do not know why it dose not work.

[Jan]

I've posted a method with a 4D array already. Using arrays instead of loops is called "vectorizing". This is sometimes very efficient in Matlab, but in this case the creation of large intermediate arrays slows doen the processing. If your supervisor asks you for a fully vectorized version (see my code above), you will learn, that a loop is more efficient in this case.

If you want to find out, why your approach does not work, post the code and a copy of the complete error message.

[Parisa Salemi]

Hi Jan. Thanks again. The problem is p1 is 3d matrix. While it must be 2D.

clc;
close all;
clear;
a=1;
c=1.5; % mm/us
T=20; dt=0.8; % us
fc=1;BW=3; % MHz
R=10;
xf=20;zf=20;
td=0;
dx=0.3; dz=0.3;pitch=0.3;
% x0=10; z0=0;
f=@(x) exp(-x.^2*BW^2);
g=@(x) sin(2*pi*fc*x);
H=@(x) f(x).*g(x);
Nz=120;
Nx=120;
x=[1:Nx]*dx;
z=[1:Nz]*dz;
xl=[1:Nx]*dx;
[z3,x3,xxl] = ndgrid(z,x,xl);
d3=sqrt((x3-xxl).^2+(z3).^2);
df=sqrt((xxl-xf).^2+(zf).^2);
td=(R-df)./c;
tt_hat=(d3./c)+td;
dim = [1, size(tt_hat)];
p   = H((0:dt:T).' - reshape(tt_hat, dim)) ./ reshape(sqrt(d3), dim);
p1  = sum(p, 4).^2;
imagesc(p1);

[Jan]

Hi. This code produces a 3D array. Of course it does, because the original code produces a set of 2D arrays. Concatenating a bunch of 2D matrices must create a 3D array.

Why do you want to get a 2D matrix? What should happen with the additional information?