C++

Forum

 
How to implement MATLAB's colon operator in C++?

JamieAl (236)
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.

MATLAB code example:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 
L = [8 8];  
Lx = L(1);
Ly = L(2);
n = [8 8];
Nx = n(1);
Ny = n(2);
xi_x =  2*pi/Lx;
xi  = (0:Nx-1)/Nx*2*pi;
x =  xi/xi_x;

yi_y =  2*pi/Ly;
yi  = (0:Ny-1)/Ny*2*pi;
y =  yi/yi_y;
% Make spatial mesh
[X,Y]   = meshgrid(x,y); 
A = 2 * pi / Lx;
B = 2*pi / Ly;

uBar = sin(A*X) .* cos(B*Y);
uk =  fft2(uBar);
n = ones(size(uBar));
nk = fft2(n);
%take convolution
productK = convolve2D(uk,nk);

where convolve2D function performs FFT product with zero padding and it looks like:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 
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
 
static const int nx = 8, ny = 8;
static const int 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!
kigar64551 (814)
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

Doesn't this boil down to [a:b c:d], which creates a vector (array) of:
[a, a+1, a+2, ..., b, c, c+1, c+2, ..., d]

(assuming that a <= b and that c <= d)

If so, it should be fairly simply to fill the array, using two loops, one that runs from a to b, and one that runs from c to d.

These two loops nested inside an "outer" loop that runs i from 1 to 2.

Of course, in each iteration of the "outer" loop, the bounds a, b, c and d need to be determined, based on current i.
Last edited on
deleted account xyzzy (5768)
If you want zero padded output look to either C's printf or C++20's std::format.

C's printf:
1
2
3
4
5
6
 
#include <stdio.h>

int main( )
{
   printf("%06d\n", 42);
}
000042

std::format:
1
2
3
4
5
6
7
 
import <format>;
import <iostream>;

int main( )
{
   std::cout << std::format("{:06}", 42) << '\n';
}

std::format has a lot of format specifiers available.

https://en.cppreference.com/w/cpp/utility/format/formatter#Standard_format_specification

If your compiler can't handle C++20, or you don't want to use C++20, there is the 3rd party C++ {fmt} library.

https://fmt.dev/latest/index.html

The <format> library is based on the {fmt} library.

std::cout can be manipulated to do padding using std::fill and std::setw.

https://en.cppreference.com/w/cpp/io/manip/setfill
https://en.cppreference.com/w/cpp/io/manip/setw

Both reside in <iomanip>, make sure you include that header.
Topic archived. No new replies allowed.