How to take the fft of a sum of exponentials?

A is a vector, and I want to take the sum of the exponential of -t/A_i over each element of A, then take the FFT. I wasn't sure what range to take for t, though it should go over several orders of magnitude. Would there be a better way to do this than just generating an arbitrary vector? Could I base it off the values of A? Could the FFT be done on a log-spaced vector? And finally, how would I figure out the frequency vector of the resulting transform? This is my attempt below.

t = linspace(0.01,100,100000);
f_t = sum(exp(-t./A));

n = 2^14;
pts = n/2; % Nyquist
s = t(2)-t(1); % linear spacing in time
Fs = (1/s)/n; % frequency to sample
Frequency = Fs*(0:pts-1);

f_w = fft(f_t,n); % fft

[Paul]

Can there be more clarification on the definition of f(t). Suppose we have A1 and A2.

The function f(t) = exp(-t/A1) + exp(-t/A2), -inf < t < inf, doesn't have a Fourier transform, at least not without some additional constraints. For example:

Is f(t) = 0 for t < 0?

Are A1 and A2 both positive?

These are basically the same point that @Matt J had about defining how f(t) behaves for large postiive/negative values of t.

Why is it desired to use FFT? If the signal has a Fourier transform, I'm going to guess it will be straightforward to derive analytically.

Answer 1

[Star Strider]

I’m not certain what you want to do. The sum function produces a scalar result if ‘t’ and ‘A’ are the same size, so the Fourier transform is a straight line.

However if ‘t’ is a row vector and ‘A’ is a column vector (making the sum argument a matrix), the ‘f_t’ result is a vector and a bit more interesting:

N = 1000;
t = linspace(0.01,100,N);
L = numel(t);
A = rand(size(t));                                          % Create 'A'
f_t = sum(exp(-t./A(:)));

% n = 2^14;
n = 2^nextpow2(L);
pts = n/2; % Nyquist
s = t(2)-t(1); % linear spacing in time
Fs = (1/s); % frequency to sample
Frequency = (Fs/2)*(0:pts-1);
ixv = 1:numel(Frequency);

f_w = fft(f_t,n); % fft

figure
plot(Frequency, abs(f_w(ixv))*2)
grid
xlabel('Frequency')
ylabel('Amplitude')
title('Linear Axes Scaling')

figure
semilogx(Frequency, mag2db(abs(f_w(ixv))*2))
grid
xlabel('Frequency')
ylabel('Amplitude (dB)')
title('Logarithmic Axes Scaling')

I made a few changes to your original code to make this work.

.

[L'O.G.]

Thanks. @Star Strider To make sure I am clear, I am trying to take the FFT of f(t) = ∑t/t_A over all A, so I think what you suggested is along those lines. Could you please tell me why you take abs(f_w(ixv))*2 ? And how to increase the range of the frequency so that I cover more decades in frequency? Thank you!

[Star Strider]

‘Could you please tell me why you take abs(f_w(ixv))*2 ?’

Sure. The fft function creates a two-sided Fourier transform, dividing the original signal energy at each frequency by 2, half allotted to the negative frequencies (after using fftshift that I do not do here because it is not necessary) and half to the positve frequencies. Multiplying by 2 in a one-sided Fourier transform corrects the signal amplitude from that allocation.

My pleasure!

Answer 2

[Matt J]

If your exponentials span too many orders of magnitude, you will probably have numerical problems, but otherwise, this is how I would set things up,

%% sampling parameters
n = 2^14;
dt=100/100000; %time sampling interval
df=1/dt/n;     %frequency sampling interval

%% set up axes
nAxis=(0:n-1)-ceil((n-1)/2);
t=nAxis*dt; %time axis
freq=nAxis*df;  %frequency axis


%% build/transform the signals
f_t = sum(exp(-t./A(:)),1);
f_w = fftshift( fft( ifftshift(f_t))  ); % fft

%% plot
plot(freq,abs(f_w))

[L'O.G.]

Thanks! Two quick things: could you tell me how you recommend changing the range of t? Also, why does your fft have the form that it does (ifftshift, then fft, then fftshift) rather than just take the fft directly?

[Matt J]

could you tell me how you recommend changing the range of t?

Did I recommend changing the range of t?

You want to choose it large enough that the exponentials have nearly decayed to zero both for positive and negative t. Come to think of it, a sum of exponentials can't decay in both directions at once, so you'll need to clarify how the signal is intended to behave for both large negative and large positive t.

why does your fft have the form that it does (ifftshift, then fft, then fftshift) rather than just take the fft directly?

The signal, as constructed in the code, has the time origin at the center of the array f_t, whereas the fft() expects the time origin to be at f_t(1). The fftshift and ifftshift just shift the origin as needed.