Row major and Column major Explained, using Python against Matlab

In both cases, a 2d array of 2 rows, 4 columns is created.
Then, it is reshaped to 1 row, 8 columns.
This clearly demonstrates how Matlab stores the array data column-wise, and Python stores the array data row-wise after it is flattened out to 1d.

Matlab (Column Major)

Screen Shot 2017-06-04 at 10.51.20 AM

>> a = = [1 2 3 4 ; 5 6 7 8]

a = 
 1 2 3 4
 5 6 7 8

a = a(:)

a =

1 
5 
2 
6 
3 
7 
4 
8

Python (Row-major)

Screen Shot 2017-06-04 at 10.51.13 AM

>> import numpy as np
>> a = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])
>> a
Out[]: 
array([[1, 2, 3, 4], 
        5, 6, 7, 8]])

>> a = a.tostring() # shows the memory arrangement, in a string

Out[] : b'\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00\x05\x00\x00\x00\x06\x00\x00\x00\x07\x00\x00\x00\x08\x00\x00\x00'

If you look at the string printed, it shows that the elements are arranged in the fashion 1, 2, 3, 4, 5, 6, 7, 8.

Why is this important to know, you ask?

For optimization, it is important. In Python, C/C++, the elements will be laid out contiguously row-wise in memory. So you should strive to process the elements row-wise. The entire row could be copied to cache and worked on from there by the OS. If you processed the elements column wise, contrary to how it is laid out in memory, you would incur expensive context switching as the OS copies in the elements you require into cache for processing for every column-wise element you process!

 

More Efficient ifftshift / fftshift in C++

matlablogo

Previously I touched upon ifftshift and fftshift in this post. https://kerpanic.wordpress.com/2016/01/15/matlab-circshift-equivalent-in-c-c/ . I’ve come to realize that there are better ways to implement fftshift and iftshift using simple memory swapping. Note that you have to pre-allocate the appropriate amount of memory for T*out before calling this function.

//-- Does 1D fftshift 
template<typename T>
inline void fftshift1D(T *in, T *out, int ydim)
{
 int pivot = (ydim % 2 == 0) ? (ydim / 2) : ((ydim - 1) / 2);
 int rightHalf = ydim-pivot;
 int leftHalf = pivot;
 memcpy(out, in+(pivot), sizeof(T)*rightHalf);
 memcpy(out+rightHalf, in, sizeof(T)*leftHalf);
}

//-- Does 1D ifftshift
//-- Note: T* out must already by memory allocated!!
template<typename T>
inline void ifftshift1D(T *in, T *out, int ydim)
{
 int pivot = (ydim % 2 == 0) ? (ydim / 2) : ((ydim + 1) / 2);

 int rightHalf = ydim-pivot;
 int leftHalf = pivot;
 memcpy(out, in+(pivot), sizeof(T)*rightHalf);
 memcpy(out+rightHalf, in, sizeof(T)*leftHalf);
}

 

Hope this helps! It is good as  a Matlab C++ equivalent.

If you want to do 2D you will have to fftshift first, transpose the matrix, then fftshift it again. Then transpose back to the original matrix form.

Also, a better circshift would be to use C++’s in-built std::rotate

template<typename T>
inline void circshift1D_IP(T *in, int ydim, int yshift)
{
 if (yshift == 0)
 return;

 if (yshift > 0) // shift right
 std::rotate(&in[0], &in[ydim - yshift - 1], &in[ydim - 1]);
 else if (yshift < 0) // shift left
 {
 yshift = abs(yshift);
 std::rotate(&in[0], &in[yshift], &in[ydim - 1]);
 }

 return;
}

Precision : Of Singles and Doubles

quote-donald-rumsfeld-be-precise-a-lack-of-precision-is-2916

Typically, a single precision digit (aka float32) has a decimal digits precision from 6-9 digits. If you type in:

eps('single')

In Matlab, you’ll find that it gives 1.1920929e-07. This is 2^-23.
Put simply, this is the smallest value required to to effect the next largest single-precision number, when the number is 1.

So for 980e6 for example, the smallest value required to represent the next largest single-precision number will be

980e6 * eps('single') == 1.1682510e+02.

Basically, if you have a single-precision value 980e6, the smallest confirmed granularity is ~116.8!  So when you are doing calculations involving large numbers, it will be good to keep this in mind. When you have something like this:

single( 980e6) + single(k(1))

single(k(1)) needs to be > ~116.8 in order to effect change upon single(980e6).
Basically I’ve said the same thing several different ways.

eps(‘double’) is 2^-52

which is much smaller.

For double-precision, it would be 980e6 * eps(‘double’) == 2.176037128265307e-07. Pretty small. Not something to fret about.