"Frequency" shift in discrete FFT in MATLAB

Question

(Disclaimer: I thought about posting this on math.statsexchange, but found similar questions there that were moved to SO, so here I am)

The context:

I'm using fft/ifft to determine probability distributions for sums of random variables. So e.g. I'm having two uniform probability distributions - in the simplest case two uniform distributions on the interval [0,1].

So to get the probability distribution for the sum of two random variables sampled from these two distributions, one can calculate the product of the fourier-transformed of each probabilty density. Doing the inverse fft on this product, you get back the probability density for the sum.

An example:

function usumdist_example()

    x = linspace(-1, 2, 1e5);
    dx = diff(x(1:2));
    NFFT = 2^nextpow2(numel(x));

    % take two uniform distributions on [0,0.5]
    intervals = [0, 0.5;
                 0, 0.5];
    figure();
    hold all;
    for i=1:size(intervals,1)
        % construct the prob. dens. function
        P_x = x >= intervals(i,1) & x <= intervals(i,2);
        plot(x, P_x);

        % for each pdf, get the characteristic function fft(pdf,NFFT)
        % and form the product of all char. functions in Y
        if i==1
            Y = fft(P_x,NFFT) / NFFT;
        else
            Y = Y .* fft(P_x,NFFT) / NFFT;
        end
    end
    y = ifft(Y, NFFT);

    x_plot = x(1) + (0:dx:(NFFT-1)*dx);
    plot(x_plot, y / max(y), '.');

end

My issue is, the shape of the resulting prob. dens. function is perfect. However, the x-axis does not fit to the x I create in the beginning, but is shifted. In the example, the peak is at 1.5, while it should be 0.5.

The shift changes if I e.g. add a third random variable or if I modify the range of x. But I can't get figure how.

I'm afraid it might have to do with the fact that I'm having negative x values, while fourier transforms usually work in a time/frequency domain, where frequencies < 0 don't make sense.

I'm aware I could find e.g. the peak and shift it to its proper place, but seems nasty and error prone...

Glad about any ideas!

Answer 1

The problem is that your x origin is -1, not 0. You expect the center of the triangular pdf to be at .5, because that's twice the value of the center of the uniform pdf. However, the correct reasoning is: the center of the uniform pdf is 1.25 above your minimum x, and you get the center of the triangle at 2*1.25 = 2.5 above the minimum x (that is, at 1.5).

In other words: although your original x axis is (-1, 2), the convolution (or the FFT) behave as if it were (0, 3). In fact, the FFT knows nothing about your x axis; it only uses the y samples. Since your uniform is zero for the first samples, that zero interval of width 1 is amplified to twice its width when you do the convolution (or the FFT). I suggest drawing the convolution on paper to see this (draw original signal, reflected signal about y axis, displace the latter and see when both begin to overlap). So you need a correction in the x_plot line to compensate for this increased width of the zero interval: use

x_plot = 2*x(1) + (0:dx:(NFFT-1)*dx);

and then plot(x_plot, y / max(y), '.') will give the correct graph: