How to vectorize loop of fft ?

7 visualizaciones (últimos 30 días)
Ole
Ole el 28 de Feb. de 2015
Comentada: Geoff Hayes el 1 de Mzo. de 2015
How to vectorize a for loop of fft ?
clear all
W = ones(1,100)*1.5;
y(1:200,1:100) = 0.3;
y(1,:) = ones(1,100)*1.2;
for j=1:200;
y(j+1,:) = ifft(fft(y(j,:)).*W);
end

Respuestas (1)

Geoff Hayes
Geoff Hayes el 28 de Feb. de 2015
Ole - is this something that you need to vectorize or just simplified? If the above is a valid example with all elements of W being identical, then the line
y(j+1,:) = (ifft((((fft((y(j,:)))).*W))));
can be simplified to
y(j+1,:) = ifft(fft((y(j,:))).*W);
= ifft(1.5*fft((y(j,:))));
= 1.5*y(j,:);
due to the linearity of the inverse Fourier transform. And so your for loop becomes
for j=1:200;
y(j+1,:) = 1.5*y(j,:);
end
The above implies that the nth column of y (for n>1) is simply
(1.5)^(n-1)*y(1,:)
in which case you can replace your code with just
y = 1.5*ones(201,100);
y(1,:) = 1.2;
y = cumprod(z);
Try the above and see what happens!
  2 comentarios
Ole
Ole el 28 de Feb. de 2015
Editada: Ole el 28 de Feb. de 2015
Thank you. W is actually a function in fourier space, I just tried to make it simple. This the structure of a fft propagation, with W being the transfer function of free space exp(i*kz*dz).
y(2:201,:) = ifft(fft(y(1:200,:)).*W); instead of the for loop gives errors
Geoff Hayes
Geoff Hayes el 1 de Mzo. de 2015
Ole - what is the error? Is it
Error using .*
Matrix dimensions must agree.
because of the line
ifft(fft(y(1:200,:)).*W);
where W is just a 1x100 array and fft(y(1:200,:)) is a 200x200 array. You can only perform an element-wise multiplication on arrays/matrices of the same dimension.
A question that I should have asked sooner - why are you trying to vectorize this problem? Performance reasons? Because the recurrence relation (k+1 depends on k) in the above for loop may provide the best solution.

Iniciar sesión para comentar.

Categorías

Más información sobre Fourier Analysis and Filtering en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by