2
$\begingroup$

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.

$\endgroup$
8
  • 2
    $\begingroup$ Please don’t cross-post. stackoverflow.com/q/77508781/7328782 $\endgroup$ Nov 19, 2023 at 0:58
  • 2
    $\begingroup$ I’m voting to close this question because cross posted here: stackoverflow.com/questions/77508781/… $\endgroup$
    – Jdip
    Nov 19, 2023 at 5:32
  • 3
    $\begingroup$ 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? $\endgroup$ Nov 19, 2023 at 15:21
  • $\begingroup$ @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 $\endgroup$
    – Engineer
    Nov 23, 2023 at 22:56
  • $\begingroup$ @Engineer, It is a great decision by MATLAB developers to make conv() based on the spatial operation and not in the frequency domain. See dsp.stackexchange.com/questions/52760. In many cases the spatial implementation is faster. I prefer the functions to be explicit and not implicit. $\endgroup$
    – Royi
    Nov 26, 2023 at 10:34

2 Answers 2

2
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ This will result in full convolution. The OP is after same. $\endgroup$
    – Royi
    Nov 26, 2023 at 10:31
  • $\begingroup$ I answered as a "hint" because this smelt like a homework problem $\endgroup$
    – Engineer
    Nov 26, 2023 at 11:38
1
$\begingroup$

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).

$\endgroup$
1
  • $\begingroup$ Any feedback form the one who -1? $\endgroup$
    – Royi
    Nov 26, 2023 at 10:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.