How to vectorize loop of fft ?

4 Ansichten (letzte 30 Tage)
Ole
Ole am 28 Feb. 2015
Kommentiert: Geoff Hayes am 1 Mär. 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

Antworten (1)

Geoff Hayes
Geoff Hayes am 28 Feb. 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 Kommentare
Ole
Ole am 28 Feb. 2015
Bearbeitet: Ole am 28 Feb. 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 am 1 Mär. 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.

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Fourier Analysis and Filtering finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by