Spectral Factorization by Radon transform

Hello,

I am interested in doing a spectral factorization of a 2D signal. To get that I am using the radom transform (radon.m and iradon.m matlab functions) to use a 1D spectral factorization instead of a 2D. I have read in the literature that you can do the radon transform, make the 1D spectral factorization to all the projections and do the iradon (tha radon transform back) to get the same result as if you had made a 2D spectral factorization.

The problem is that I am not getting very good results and I do not know if I am doing something wrong (the spectral factorization 1D function is totally checked).

Thanks in advance

[Sean de Wolski]

Can you post your code?

Answer 1

[Matt J]

I'm assuming that by "spectral factorization" you mean "Fourier Transform".

Then, no. IRADON wouldn't be involved. The 1D spectra form the radial cross-sections of the 2D spectra through the 2D frequency origin. In other words, they give you samples of the 2D spectra on a polar coordinate grid. You would use interp2 or something to interpolate back to Cartesian coordinates.

The result will be influenced by interpolation artefacts as well as artefacts in the output of RADON, which uses a rather crude algorithm, unfortunately.

Why exactly are you interested in doing it this way? It's much more efficient and accurate to do a 2D Fourier transform directly.

Answer 2

[Peter]

Hello,

Thanks for your answers.

No, I actually mean Spectral Factorization. If you have an autocorrelation signal you can compute the minimum-phase impulse response which satisfies that given auto-correlation.

If you have 1D signal, this process is quite simple. If you have a multidimensional signal this process is not simple at all. Then, you can try to do a 2D spectral factorization (wich solution is not guaranteed) or you can try to reduce the 2D signal to an equivalent 1D signal and solve it with 1D processes.

I have seen two methods to "reduce" the dimensionality. First of them is the Helical boundary and the second is the Radon projection. (I think the helical boundary is a special case of the radon projection)

A simple code to do that would be:

 theta=1:180;
 [r_theta,xp] = radon(r_x,theta,M);    
 for ii=1:length(theta)
        h_theta(:,ii)=SpectralFactorization(r_theta(:,ii));
 end
 SF=iradon(real(h_theta),theta,'linear','None')

being theta=points where I am doing the projection, r_x= the 2D autocorrelation signal, r_theta= All the radon projections, h_theta the spectral factorizations of the radon projections, SF the spectral factorization of a 2D signal obtained with 1D processes.

...but the solution is not as good as I know that it could be.

Any ideas,

Thanks.

Answer 3

[Peter]

Any thoughts?

[Matt J]

FYI, the "Answers" in these threads are not meant as a means of dialogue. You should respond to other peoples' Answers using Comments. If you have follow-up to your original Question, you can post that as a Comment to that Question, like Sean has done.

Aside from that, neither the problem you are attacking, nor the technique used to solve it look like common knowledge. It might help if you post some links to relevant sources. In particular an impulse response that "satisfes a given auto-correlation" seems like weird lingo to me.

It could be that the artifacts you are seeing are purely due to under-sampling. Maybe if you applied radon over more angles, it would improve the result. As I also mentioned earlier, RADON uses a rather crude algorithm. The Siddon algorithm, e.g., as used here

http://www.mathworks.com/matlabcentral/fileexchange/24479-pet-reconstruction-system-matrix

might give better results.

[Peter]

Ok, I did not know that policy(about comments and answers).

It is not about under-sampling. I have tried with so many angles and the result is pretty the same.

"satisfes a given auto-correlation" means that normally you have a signal and then you can obtain the autocorrelation. Here is different, you have an autocorrelation and you would like to know which signal correspond to that autocorrelation sequence.

I know it is not common knowledge but such I am using Matlab functions and there is some people who usually use them, I think this web was my better choice.

Thanks anyway