I have the following MATLAB code that uses zero padding to deal with aliasing errors that I would like to rewrite in my C/C++ code. The issue is with the zero padding technique is easily implemented with MATLAB's operator (:) unlike in the C/C++ code.
function hk = convolve2D(fk, gk)
% Takes convolution of two 2D matrices in Fourier space
% Pad f and g with zeros.
N = size(fk);
% Pad with zeros to get a total of 3/2 N to get optimal computational efficiency
M = 3 * N / 2;
% Make a new set of matrices of only zeros.
fk_pad = zeros(M(1), M(2));
gk_pad = zeros(M(1), M(2));
% Get indices of the parts of the matrix we want to populate
for i = 1:2
ind_pad{i} = [1:N(i)/2 M(i)-N(i)/2+1:M(i)];
end
% Populate zero matrices with fk and gk
fk_pad(ind_pad{1},ind_pad{2}) = fk;
gk_pad(ind_pad{1},ind_pad{2}) = gk;
% Take convolution as before. Note we need to renormalize to the originally sized FFT of N
hk_pad = (prod(M) / prod(N)) * fft2(real(ifft2(fk_pad) .* ifft2(gk_pad)));
% Choose only the indices we care about for hk
hk = hk_pad(ind_pad{1}, ind_pad{2});
% Get rid of additionally padded zeros for the highest order modes
hk(N(1)/2+1,:) = 0;
hk(:,N(2)/2+1) = 0;
end
The particular line of code I am struggling to implement in C++ is:
1 2 3
for i = 1:2
ind_pad{i} = [1:N(i)/2 M(i)-N(i)/2+1:M(i)];
end
Is there a way for me to do this in C++ using a for loop only and without the use of any special matrix class. What I currently have in C++ (which is not much, basically initializing my arrays/matrices):
1 2 3 4 5 6 7
staticconstint nx = 8, ny = 8;
staticconstint Mx = 3*nx/2, My= 3*ny/2;
double N[2] = {nx,ny};
// Make a new set of matrices of only zeros.
fftw_complex *fk_pad;
fk_pad = (fftw_complex*) fftw_malloc((Mx*My)* sizeof(fftw_complex));
memset(fk_pad, 42, (Mx*My) * sizeof(double)); //memset 42 initializes array to 1.42603e-105
I am using fftw_malloc here on purpose since I am heavily relying in FFTW3 library. Thanks!