I have 46 MODIS NBAR images of the same area (time-series). How can I get seven (1 for each band) multilayer (46 layer) images, for the further FastICA analysis? According to the results of the paper, it should all result in matrix A of 46*10 in size and 10 independent source images. Input to FastICA : NDVI or SR.
you can do this with gdal, more specifically gdalbuildvrt.
first you can build vrt for each band and each image (-b to select the band number), then you can assemble those bands together (-separate to stack images together). Finally, you can use gdal_translate to export the vrt to a new format (if needed). Further details on how to use MODIS hdf can be found on this site, e.g. extract scientific layers from MODIS HDF Dataeset using Python GDAL