matlabimage-processingfftwindowed

MATLAB windowed FFT for a 2D matrix data (image)


Data: Say I have a 2000 rows by 500 column matrix (image)

What I need: Compute the FFT of 64 rows by 10 column chunks of above data. In other words, I want to compute the FFT of 64X10 window that is run across the entire data matrix. The FFT result is used to compute a scalar value (say peak amplitude frequency) which is used to create a new "FFT value" image.

Now, I need the final FFT image to be the same size as the original data (2000 X 500).

What is the fastest way to accomplish this in MATLAB? I am currently using for loops which is relatively slow. Also I use interpolation to size up the final image to the original data size.


Solution

  • As @EitanT pointed out, you can use blockproc for batch block processing of an image J. However you should define your function handle as

    fun = @(block_struct) fft2(block_struct.data);
    B = blockproc(J, [64 10], fun);
    

    For a [2000 x 500] matrix this will give you a [2000 x 500] output of complex Fourier values, evaluated at sub-sampled pixel locations with a local support (size of the input to FFT) of [64 x 10]. Now, to replace those values with a single, e.g. with the peak log-magnitude, you can further specify

    fun = @(block_struct) max(max(log(abs(fft2(block_struct.data)))));
    B = blockproc(J, [64 10], fun);
    

    The output then is a [2000/64 x 500/10] output of block-patch values, which you can resize by nearest-neighbor interpolation (or something else for smoother versions) to the desired [2000 x 500] original size

    C = imresize(B, [2000 500], 'nearest');
    

    I can include a real image example if it will further help.

    Update: To get overlapping blocks you can use the 'Bordersize' option of blockproc by setting the overlap [V H] such that the final windows of size [M + 2*V, N + 2*H] will still be [64, 10] in size. Example:

    fun = @(block_struct) log(abs(fft2(block_struct.data)));
    V = 16; H = 3; % overlap values
    overlap = [V H]; 
    M = 32; N = 4; % non-overlapping values
    B1 = blockproc(J, [M N], fun, 'BorderSize', overlap); % final windows are 64 x 10
    

    However, this will work with keeping the full Fourier response, not the single-value version with max(max()) above.

    See also this post for filtering using blockproc:Dealing with “Really Big” Images: Block Processing.