My goal (to spare space and computing time) is then to transform the sensor-space TF map using directly the imaging kernel A. Since I use a current density estimation (MNE) and the Fourier transform is linear this (I think) is possible.
Yes, this works. This solution is offered in the process options when computing time-frequency maps on sources ("Optimize storage: Yes, save morlet(sensors) and inverse kernel").
https://neuroimage.usc.edu/brainstorm/Tutorials/TimeFrequency#Full_cortical_maps
The next step is to select only a few scouts from the source space projected complex TF map. The problem is that the atlas I'm currently using is defined only on the default anatomy.
Where did you get your subject anatomy from? If it was processed with FreeSurfer or CAT, you should have the same atlases available in your subject and in the MNI template. If it was not processed with FreeSurfer or CAT, you don't have any reliable way to get a correspondence between your subject anatomy and your template.
https://neuroimage.usc.edu/brainstorm/Tutorials/CoregisterSubjects
I have to project (in some way) the source TF map to the default anatomy
You could also project the scouts you need on the individual anatomy. This is a lot lighter solution.
If you don't need full source time-frequency maps on the MNI template for group analysis, don't compute them at all.
Where I can find, or how can I compute, the interpolation matrix?
The interpolation matrix is computed in function tess_interp_tess2tess.m. This is not easy to use manually... See how this is done from bst_project_sources.m and bst_project_scouts.m