How can I evaluate characteristic functions in MatLab?

I want to numerically evaluate characteristic functions (CF) of PDFs (definition: https://en.wikipedia.org/wiki/Characteristic_function_(probability_theory) )

For example, the CF of the Normal distribution is

where t is the CF variable, μ is the distribution mean, is its variance and i is the imaginary unit. If I numerically construct a normal distribution as below, how could I numerically estimate its charactertic function?

avg = 1; % average
s   = 1; % standard deviation

dx = 0.01;
x = (avg - 5*s):dx:(avg+5*s);

gaussian = 1/sqrt(2*pi*s^2) * exp ( - 0.5 *( x - avg ).^2 / s^2 ); 

EDIT: I have found an answer that works if one has access to the data generated by the PDF. In the case of a Gaussian distribution with average avg and standard deviation s, one could estimate the CF by the Empirical Characteristic Function (ECF, detailed in Feuerverger & Mureika 1977):

In code,

N = 1e4; % Number of data points
y = m + s*randn(1,N); % N gaussian random numbers with average avg and std s

u = 0.1:0.05:3; % Argument of the characteristic function

ECF = zeros(1,length(u)); % Emperical characteristic function

for j = 1:length(u)
    
    ECF(j) = 1/N * sum( exp ( 1i * u(j) * y ) );
    
end

The answer provided by Paul below yields equivalent results.

Answer 1

[Paul]

The CF is the Fourier transform of the pdf. Approximate the pdf as a sum of line segments. Take the sum of the Fourier transforms of each segment.

syms g(x) a b c m w
assume([x a b c m w],'real');
g(x,a,b,c,m) = rectangularPulse(a,b,x)*(m*(x-a)+c);    % equation of line segment
G(w,a,b,c,m) = simplify(conj(fourier(g(x,a,b,c,m))));  % Fourier transform of line segment, conj required when using the default FourierParameters

Gfunc = matlabFunction(G)

%{
Gfunc =

  function_handle with value:

    @(w,a,b,c,m)-1.0./w.^2.*(m.*exp(a.*w.*1i)-m.*exp(b.*w.*1i)-c.*w.*exp(a.*w.*1i).*1i+c.*w.*exp(b.*w.*1i).*1i-a.*m.*w.*exp(b.*w.*1i).*1i+b.*m.*w.*exp(b.*w.*1i).*1i)
%}

% example with standard normal distribution

% distribution parameters
mu = 1;
sigma = 1;

% the pdf
x = -4:.1:4;
gpdf = normpdf(x,mu,sigma);

% generate the CF
w = -4:.01:4;
gcfunc = 0*w;
for ii = 1:numel(x)-1
    
    gcfunc = gcfunc + Gfunc(w,x(ii),x(ii+1),gpdf(ii),(gpdf(ii+1)-gpdf(ii))/(x(ii+1)-x(ii)));
    
end
    
% true CF for comparison
cftrue = (exp(1i*mu*w - sigma^2*w.^2/2));

% compare
subplot(211);plot(w,real(gcfunc),w(1:10:end),real(cftrue(1:10:end)),'o'),grid
subplot(212);plot(w,imag(gcfunc),w(1:10:end),imag(cftrue(1:10:end)),'o'),grid

[Martim Zurita]

Paul, thanks so much for your answer. Not only it solved my problem, it also deepened my knowledge of programming and Fourier analysis.

[Paul]

Glad you found it helpful.

I don't know how much the tails of a distribuiton, should they exist, contribute to the CF, or how close the spacing needs to be. Maybe there's a way to check some properties of the computed CF to determine if the answer is reasonable.

Answer 2

[Jeff Miller]

It seems like your question contains its own answer:

mu = 1;
sigma = 1;
t = -2.5:0.01:2.5;
cf = exp(i*mu*t - sigma^2*t.^2/2);
plot(t,cf)

Or maybe I'm missing something--do you want to be able to do this for any pdf, even when you don't have a handy expression for the cf?

[Martim Zurita]

Sorry Jeff, I may not have been clear on my question. Hopefully Paul was still able to grasp what I meant. Nevertheless, thanks for your attention.