how to raise the fourier transformed probability density function to a fractional power?

Question

By hypothesis my measured probability density functions (PDF) result from n convolutions of an elementary distribution (E).

I have two distributions the first (F) of which is supposed to have undergone more convolutions (m_1) than the second (G) (m_2 convolutions).

In fourier space:

F' = E'^m_1
G' = E'^m_2

As the two PDFs are constituted from the same elementary distribution, I should be able to be able to calculate the PDF of G from F

G' = F'^{m_1/m_2}

Taking the IFFT i should have a distribution that overlaps well with G.

A naive way of doing this would to be simply to calculate the FT of F and raise it to the power 1/integer and testing a range of integers.

My question are there any tricks for raising the Fourier transformed PDF to a fractional power. I have done so but the IFFT gives a distribution far from that which is expected. And strange aliasing errors.

I've included a padded vector as one might do if they were to do a convolution of two PDFS.

My normalization is based on the fact that the k=0 [ProbF(1,1)] wave vector gives the integral of the PDF which should be equal to one.

Of course, the hypothesis could be wrong, but it has all the reasons in the world to be valid.

My code

Inc  = INC1 ;                           % BINS 
null = zeros(1,length(Inc));            % PADDED PROB
Inc  = [ Inc.*-1 (Inc)  ];              % PADDED INC VECTOR
Prob = [ null heightProb1  ]       ;    % PADDED PROB VECTOR

ProbF     = (fft(Prob))         ;
ProbFnorm = ProbF./ProbF(1,1) ;         % NORMALIZED BY K=0 COMPONENT (integral of PDF =1)

m=.79                                   % POWER TO RAISE

ProbFtrans = ((ProbFnorm).^(m));        % 'DECONVOLUTION' IN FOURIER SPACE
ProbIF     = (ifft((ProbFtrans)).*(ProbF(1,1)));    % RETURN TO PROBABILITY SPACE

figure(2);
plot(Inc,ProbIF,'rs')

Thank you in advance for your help

Emerald Weapon Emerald Weapon · Accepted Answer · 2014-05-28T17:17:07

Fourier coefficients are typically complex numbers (unless your function is symmetric).

You should be very careful when you raise complex numbers to fractional powers.

For example, consider

z=1.2 + i*0.65;

then raise z to power 4

>> z^4

ans =

     -1.398293750000001e+00 + 3.174599999999999e+00i

and to power 8.

>> z^8

ans =

     -8.122859748710933e+00 - 8.878046677500002e+00i

Then try obtain z^4 as (z^8)^(1/2)

>> (z^8)^(1/2)

ans =

      1.398293750000001e+00 - 3.174600000000000e+00i

Surprise! You don't get z^4! (wrong sign)

If you avoid taking the fractional power and "rewind" z^8 by diving back by z you get back z^4 correctly:

>> (z^8)/z/z/z/z

ans =

     -1.398293750000000e+00 + 3.174599999999998e+00i

The reason is in the definition of fractional powers in the complex plane. Fractional powers are multi-valued functions, which are made single-valued by introducing a branch-cut in the complex plane. The nth-root z^(1/n) has n possible values, and
matlab singles out one of these by interpreting the complex function z^(1/n) as its so-called principal branch. The main implication is that in the world of complex numbers ^1/n does not always invert ^n.

If this doesn't make any sense to you, you should probably review some basic complex analysis, but the bottom line is that fractional powers of complex numbers are tricky animals. Wherever possible you should try to work around fractional powers by using division (as show above).

I am not sure this will fix all of your problems, but from your description it looks like this is certainly one problem you have.

how to raise the fourier transformed probability density function to a fractional power?

1 Answers