Matlab: Setting complex-number array values is slow

Question

I'm a Matlab Newb. I'm currently experimenting with Fourier Transformation and trying to implement a frequency filter (deleting all frequencies but those between k_min and k_max).

In order to realize that, I'm deleting the respective pixels in the Fourier transformed image. I'm using the following code:

% Example values
kmin = 0;
kmax = 300;

for i = 1:w
    for j = 1:h
        if norm([w/2, h/2] - [i,j]) < kmin || norm([w/2, h/2] - [i,j]) > kmax
            Fs(j,i) = 0.0;
        end
    end
end

My image is about 1000x600. Fs is therefore a 1000x600 array of COMPLEX NUMBERS. Now my problem: Why is this so slow? If I set very few pixels to zero (e.g. kmin = 10, kmax = infinite) then the code runs quickly but if I have to set almost all pixels to 0.0 (e.g. kmin = 0, kmax = 10) it takes an incredible amount of time to complete.

All I do is set some array entries to zero (worst case less than 1,000,000 entries, maybe factor two because they are complex numbers). Why does this take up to minutes? :)

Answer 1

Loops tend to be slow in Matlab (although that's improving in recent versions). You better vectorize, which means working with all values "at the same time", without loops. In this case you can do it using ndgrid and linear / logical indexing:

[ii jj] = ndgrid(1:w,1:h); %// creates all combinations
aux = sqrt((w/2-ii(:)).^2+(h/2-jj(:)).^2); %// compute norm
ind = (aux<kmin) | (aux>kmax); %// logical index
Fs(ind) = 0;

Answer 2

Instead of calculating the frequency from the index and seeing whether it falls into the passback, calculate the index corresponding to the cutoff frequency.

Once you have n_cutoff_low and n_cutoff_high calculated, setting all elements from 1:n_cutoff_low and from n_cutoff_high:end to zero is easy and very fast.

Basically, you use fftshift to change your formula from

norm([w/2, h/2] - [i,j]) < kmin

to

norm([i,j]) < kmin

then square this to

i^2 + j^2 < kmin^2

which means

j^2 < kmin^2 - i^2

or

j_cutoff = floor(sqrt(kmin^2 - i^2));
Fs(i,-j_cutoff:j_cutoff) = 0;