This function has a few issues. First you want to predict signal from the coefficients of the csd fit. You can do this by using simply taking `dot(Pmatrix, fit.shm_coeff)`. In this case `Pmatrix == dot(SH_basis * model.R)`. You can take a look at the fit method to see how we get the SH-basis and `model.R` is the spherical harmonic representation of the response function, meaning that `coef-SIG = dot(R, coef-FOD)`. You can re-compute R if you want, but it's already computed for every CSD model.

Keep in mind with the response function that it needs to be axially symmetric, this is taken care of by the `sh_to_rh` function when R is created. We've actually done something very weird with the response function, we represent it as a tensor and that's not very good. What actually happens is that an axially symmetric SH function is fit to the tensor and that's the true response function.