2

In MATLAB, I tried to convert the fft2-based multiplication using conv2.

Z = rand(100);
k = rand(100);
Fk= fft2(k);
V1 = real(ifft2(Fk.*ifftshift(Z)));
V2 = conv2(ifft2(ifftshift(Z)),k, 'same');

I want V1 and V2 to be the same, but they're not. I need to avoid using the FFT while calculating V2.

Royi
  • 19,608
  • 4
  • 197
  • 238
Ka Mirul
  • 21
  • 1
  • 2
    Please don’t cross-post. https://stackoverflow.com/q/77508781/7328782 – Cris Luengo Nov 19 '23 at 00:58
  • 2
    I’m voting to close this question because cross posted here: https://stackoverflow.com/questions/77508781/convert-fft2-version-to-conv2-matlab – Jdip Nov 19 '23 at 05:32
  • 3
    this is not surprising – real(IFFT(a .* b)) simply isn't the same as conv(a, b). The first one discards the imaginary part, the second doesn't. The first one does a cyclic convolution, the second doesn't. They are simply not the same operation. AND, big surprise, the implementation of conv2 almost surely internally uses an FFT (though of a different size than you are), because convolution is usually implemented using appropriate padding and FFTs. Which brings one to the next question: why would you want to avoid the FFT? – Marcus Müller Nov 19 '23 at 15:21
  • @MarcusMüller something interesting/weird about MATLAB is that the conv function doesn't use FFT, it actually does "flip-and-shift" the signals. Using conv with anything more than a few hundred samples becomes very slow, and so one can use the xcorr function, where FFTs are used, as a proxy for convolution – Engineer Nov 23 '23 at 22:56
  • @Engineer, It is a great decision by MATLAB developers to make conv() based on the spatial operation and not in the frequency domain. See https://dsp.stackexchange.com/questions/52760. In many cases the spatial implementation is faster. I prefer the functions to be explicit and not implicit. – Royi Nov 26 '23 at 10:34
  • @Royi yup, it is documented in a Mathworks blog here: https://blogs.mathworks.com/steve/2009/11/03/the-conv-function-and-implementation-tradeoffs/. Depends on how many samples you use! Seems to align with your previous answer too – Engineer Nov 26 '23 at 11:37
  • @Engineer, This is a very old conv(). They have improved it greatly since then. So the line higher. Nothin new in this, there is asymptotic complexity and the constant. – Royi Nov 26 '23 at 12:06
  • @Royi ok...sounds like we agree – Engineer Nov 26 '23 at 13:15

2 Answers2

2

HINT

%% 1D example
Nx = 200; % number of samples in x
Ny = 100; % number of samples in y

x = randn(Nx, 1); % random data x
y = randn(Ny, 1); % random data y


xp = padarray(x, Ny-1, 0, "post"); % pad x with 0's
yp = padarray(y, Nx-1, 0, "post"); % pad y with 0's


z_conv = conv(x, y); % direct method
z_fft = ifft(fft(xp) .* fft(yp)); % FFT based method
Engineer
  • 3,032
  • 1
  • 8
  • 16
1

You have a MATLAB Code as an answer in Replicate MATLAB's conv2() in Frequency Domain.

For discrete signals, the multiplication in frequency domain is equivalent of circulant convolution. In order to achieve linear convolution with different shapes (full, same, valid) you need to do some padding and indices arithmetic as shown in the linked answer.

Julia Implementation

As a practice, let's implement this in Julia.

using FFTW;

function Conv2DFreqDom( mI :: Matrix{T}, mH :: Matrix{T}; convMode :: ConvMode = CONV_MODE_FULL ) where {T <: AbstractFloat}


numRowsI, numColsI = size(mI);
numRowsH, numColsH = size(mH);

if (convMode == CONV_MODE_FULL)
    numRowsFft  = numRowsI + numRowsH - 1;
    numColsFft  = numColsI + numColsH - 1;
    mO = Matrix{T}(undef, (numRowsI, numColsI) .+ (numRowsH, numColsH) .- 1);
elseif (convMode == CONV_MODE_SAME)
    numRowsFft  = numRowsI + numRowsH;
    numColsFft  = numColsI + numColsH;
    mO = Matrix{T}(undef, (numRowsI, numColsI));
elseif (convMode == CONV_MODE_VALID)
    numRowsFft  = numRowsI;
    numColsFft  = numColsI;
    mO = Matrix{T}(undef, (numRowsI, numColsI) .- (numRowsH, numColsH) .+ 1);
end

mT1 = Matrix{complex(T)}(undef, numRowsFft, numColsFft);
mT2 = Matrix{complex(T)}(undef, numRowsFft, numColsFft);

Conv2DFreqDom!(mO, mI, mH, mT1, mT2; convMode = convMode);

return mO;




end


function Conv2DFreqDom!( mO :: Matrix{T}, mI :: Matrix{T}, mH :: Matrix{T}, mT1 :: Matrix{Complex{T}}, mT2 :: Matrix{Complex{T}}; convMode :: ConvMode = CONV_MODE_FULL ) where {T <: AbstractFloat}


numRowsI, numColsI = size(mI);
numRowsH, numColsH = size(mH);

if (convMode == CONV_MODE_FULL)
    numRowsFft  = numRowsI + numRowsH - 1;
    numColsFft  = numColsI + numColsH - 1;
    firstRowIdx = 1;
    firstColIdx = 1;
    lastRowIdx  = numRowsFft;
    lastColdIdx = numColsFft;
elseif (convMode == CONV_MODE_SAME)
    numRowsFft  = numRowsI + numRowsH;
    numColsFft  = numColsI + numColsH;
    firstRowIdx = ceil((numRowsH + 1) / 2);
    firstColIdx = ceil((numColsH + 1) / 2);
    lastRowIdx  = firstRowIdx + numRowsI - 1;
    lastColdIdx = firstColIdx + numColsI - 1;
elseif (convMode == CONV_MODE_VALID)
    numRowsFft  = numRowsI;
    numColsFft  = numColsI;
    firstRowIdx = numRowsH;
    firstColIdx = numColsH;
    lastRowIdx  = numRowsFft;
    lastColdIdx = numColsFft;
end

mT1[1:numRowsI, 1:numColsI] .= mI;
mT2[1:numRowsH, 1:numColsH] .= mH;
fft!(mT1);
fft!(mT2);

mT1 .*= mT2;
ifft!(mT1);

mO .= real.(@view mT1[firstRowIdx:lastRowIdx, firstRowIdx:lastRowIdx]);




end


Padding


mIPad = PadArray(mI, (kernelRadius, kernelRadius); padMode = padMode);


Generating the Reference


mORef = Conv2D(mIPad, mH; convMode =  convMode);


Using FFT


mODft = Conv2DFreqDom(mIPad, mH; convMode =  convMode);

This will yield a perfect result as:

julia> maximum(abs.(mORef - mODft))
1.1102230246251565e-15

The full Julia code is available on my StackExchange Signal Processing Q90036 GitHub Repository (Look at the SignalProcessing\Q90036 folder).

Royi
  • 19,608
  • 4
  • 197
  • 238