Transformation from spectrum into pulse

I have a chirped pulse of the form

E(t) = exp(-t^2)*cos(-15t+7t^2).

Which looks like this

If I transform this pulse into frequency domain analytically, that is manually on a piece of paper I will get the following spectrum:

E(w) = E1(w) + E2(w),

where

E1(w) = exp(-(w-15)^2/200)*exp(-7i(w-15)^2/200)

E2(w) = exp(-(w+15)^2/200)*exp(-7i(w+15)^2/200)

Now I want to check if this is really the case using ifft. I give the function to be inverse-transformed by ifft to be the spectrum with mathematical form as written above, after applying ifft to this function I hope I can recover the chirped pulse. Here's the code I write to do this:

Fs = 500;
w = -100:1/Fs:100;
Ew = exp(-(w-15).^2/200).*exp(-7*1i*(w-15).^2/200) + exp(-(w+15).^2/200).*exp(7*1i*(w+15).^2/200); % Spectrum of a chirped pulse.
nfft = 2000000;
Et = ifft(Ew,nfft); Pulse in time domain.
Et = Et(1:nfft);
Et = ifftshift(Et);

Et = real(Et);
t = (-nfft/2:(nfft/2-1))*Fs/nfft;
subplot(2,1,1)
plot(w,Ew)
subplot(2,1,2)
plot(t,Et)

And I what I got as the pulse is like this,

Which clearly not the original pulse I started with. Can somebody tell me what the problem was? I have tried to review my manual calculation and I could find no mistake in it. I omitted some prefactors but that shouldn't pose any effect on the form of the curve.

Thanks in advance for your help.

[Matthew Crema]

I think there are a number of issues here, including some confusion about when to use frequency in Hz, and when to use radian frequency in rad/sec. Second, I'm not sure why you chose the sample rate and number of points that you did. For something like this why not use the same number of samples in both the time and frequency domain? Third, the signal starts before time = 0, which (for me at least) is unusual, so I'd start by time shifting the whole signal to deal with only positive times (and be careful to remember this was done through the code). Finally, I have some doubts about the correctness of your expression for E(w), which I admit I have not played with on paper. But bear in mind that whenever you take exp(x) the quantity x should be dimensionless, which in your expression it does not seem to be.

So I wrote some code from scratch, which still doesn't get you all the way there, but may be helpful.

%%Constants

% Number of points in discrete signals [samples]
Npts = 256;
% Sampling Rate [samples/second]
Fs = 32;
% Duration of time domain signal [sec]
Dur = Npts/Fs;
% Indexing vector n for both time and freq vectors [samples]
n = (0:Npts-1);
% Time vector [sec]
time = n/Fs;
% Frequency vector freq = (0:Npts-1)./Npts*fs [Hz]
freq = n/Dur;

%%Theoretical Functions

% Compute known time domain signal (shift by tshft to make it causal)
w1 = 15;    % Parameter 1 (assumed units are rad/sec)
w2 = 7;     % Parameter 2 (assumed units are rad/sec^2)
w3 = 1;     % Parameter 3 (pulse width, assumed units are sec)
tshft = 4;  % Shift by 4 seconds so only have to deal with positive time

% Given Expression, all units check out OK
Et_orig = exp(-(time-tshft).^2./w3^2) .* ...
    cos(-w1*(time-tshft) + w2*(time-tshft).^2);

% Compute fft of Et_orig to later compare with Ew_wrong
Ew_correct = fft(Et_orig);

% Real-valued time domain signal obtained by ifft
Et_correct = ifft(Ew_correct);

%%Build Spectrum from given (but incorrect) exprssion and Compute IFFT

% Determine indices and frequency vector for single-sided
% spectra (First 1/2 of FREQ and FFT vectors)
halfpt = floor(Npts/2)+1;  % Index for Nyquist frequency
ssinds = 1:halfpt;  % Indices for single-sided spectra (f>=0)
f_half = freq(ssinds);  % s-s Freq vector (0 <= f <= fs/2)

% Compute given (possibly incorrect) spectrum of a chirped pulse
sf = 200;   % some scale factor? what are its units?
w = 2*pi*f_half;

% Given Expression. Units do not check out.
Ew_incorrect = exp(-(w-w1).^2/sf).*exp(-1j*w2*(w-w1).^2/sf) + ...
    exp(-(w+w1).^2/sf).*exp(1j*w2*(w+w1).^2/sf);

% Time shift for agreement with time domain signal
Ew_incorrect = sqrt(Npts)* Ew_incorrect.*exp(-1j*2*pi*f_half*tshft);

% Build Conjugate-symmetric sequence so IFFT returns the correct real
% valued sequence
Ew_incorrect(halfpt+1:Npts) = conj(fliplr(Ew_incorrect(2:ceil(Npts/2))));

% Real-valued (incorrect) time domain signal obtained by ifft
Et_incorrect = real(ifft(Ew_incorrect));

%%Plot
subplot(3,1,1)
plot(time, Et_orig, time, Et_correct, 'g--', time, Et_incorrect, 'r:')
ylabel('x(t)')
xlabel('Time (sec)')
subplot(3,1,2)
plot(freq, real(Ew_correct), 'g--', freq, real(Ew_incorrect), 'r:')
ylabel('Real(FFT(x))')
lg = legend('Numerical fft(x)', 'Given Fourier Transform');
set(lg, 'FontSize', 6)
subplot(3,1,3)
plot(freq, imag(Ew_correct), 'g--', freq, imag(Ew_incorrect), 'r:')
ylabel('Imag(FFT(x))')
xlabel('Frequency (Hz)')

[Matthew Crema]

Also note that it may be easier (or more difficult) to perform the analytical transform if you think of your pulse as a linear chirp with parameters f0 and f1, modulated by a Gaussian envelope with parameter Tau. In your case the values for these are given below, but for analytical treatment just treat them as constants,

% Parameters
f0 = -2.3873;  % Instantaneous Freq at beginning of chirp (t=0) [Hz]
f1 = 15.4383;   % Instanteneous Freq at end of chirp (t=Dur) [Hz]
Dur = 8; % Time for instantaneous frequency to increase from f0 to f1
k = (f1-f0)/Dur; % Chirp rate (f1-f0)/Dur
Tau = 0.5;  % Related to Width of Gaussian Envelope [sec]
tshft = 4; % Shift by 4 seconds
tnorm = time - tshft;
Phi_t = 2*pi*(f0 + k/2*tnorm).*tnorm;  % Instantaneous Phase

% Gaussian Envelope
env_t = exp(-tnorm.^2./(2*Tau)^2);

% Equivalent Expression for E(t)
E(t) = env_t .* cos(Phi_t);

Note that the wikipedia page for "Chirp" states

"In the special case where [env_t] is constrained to be a flat topped pulse ... then an analytical solution [for the Fourier transform] is possible."

This solution is given on that page. However the statement above may be hinting that a closed-form analytical solution for your function (whose envelope is a Gaussian pulse) does not exist. This would not be too surprising if you find at some point you have to integrate the Gaussian function, as the integral of a Gaussian has no closed form solution (which is why we traditionally used look up tables and now compute it numerically using an approximation).

However, my intuition suggests that you will instead end up performing not integration, but convolution with (the Fourier Transfrom of) the Gaussian, which is usually quite possible.

[Imam]

Hi Matthew, thanks for your response.

The number of samples in frequency and time domains are not necessarily the same. In fact most people did the so-called zero padding, which you might have heard of, to purposely increase the number of samples in the after-transformation domain (in my case freq domain). This is done to smooth the transformed curve, I have checked your first code and found that the last two curves have jagged sharp peaks which should have been smoothly curving. But I guess this is not the major problem in my question. I know that the argument of exponential should be dimensionless, and the ones appearing in my expression above have implicitly fulfilled that. Only that because all constants are already expressed in numbers then it's not seen. Anyway, Matlab doesn't care about whether thing has dimension or not.

Btw I'm pretty sure of my calculation since it was based on the lecture notes from a related course. Well for the sake of easiness I will replace that problem with a simpler one. Now I remove the chirp, let's say the pulse has the form

E(t) = exp(-t^2/4)*cos(4t)

This looks like

Now I do the same thing as before, namely compute manually the spectrum from using FT and I end up with

E(w) = exp(-(w-4)^2) + exp(-(w+4)^2)

I guess that was a typical integration that one working in optical telecommunication will normally encounter. The next is still the same thing, that is trying to prove the correctness of that expression using MATLAB. However I did it differently than before, now I wouldn't use MATLAB's build-in ifft function, instead I would go through the exact formula of Fourier transformation:

Since this involves integration, I employ numerical integration algorithm, trapz function, into my new code:

w = -10:1/100:10; % Frequency (rad/s)
f0 = 4; % Center frequency
Ew = exp(-(w-f0).^2) + exp(-(w+f0).^2); * Spectrum
t = -6:0.01:6; % Time (s)

for i=1:length(t)
  for j=1:length(w)
  A(i,j) = Ew(j)*exp(1i*w(j)*t(i));

    end
    Et(i) = trapz(w,A(i,:));
end

Et_R = real(Et);
subplot(2,1,1)
plot(w,Ew)
xlabel('frequency (rad/s)')
ylabel('E(w)')
subplot(2,1,2)
plot(t,Et_R)
xlabel('time (s)')
ylabel('E(t)')

Amazingly the pulse in time domain E(t) exactly represents the original pulse.

[Imam]

But if I kept using the ifft function to transform the expression of spectrum in my previous comment, what I got is

Now this makes me question what the fft and ifft function actually do. How can its result deviates from the fundamental FT equation.

[Salaheddin Hosseinzadeh]

Hi Imam,

One thing I don't understand is why you're using IFFT instad of FFT?! There are similarities and differences.

You have the time domain signal, you calculated manualy the frequency domiain by using Fourier Transform perhaps and you want to verify if you did it correct!

Alright!

1 - plot your manualy solved frequency equation in MATLAB

2 - Take FFT of the time domain signal

And compare them, that's all!

Is this not all you want to do?!

[Imam]

Hi Salaheddin, what you said is true about what I wanted to do.

Why did I use ifft instead of fft? It's because I have function in time and wanted to transform it into frequency. Well in fact the result of fft is just the mirror of ifft with respect to zero, but the fucntional form is preserved.

And I actually have tried to do what you suggested me above well ahead of time, namely compare the spectrum from manual calculation and that from MATLAB computation. The result was unsatisfying, you see the expression of spectrum I got manually in previous comments, it's a real function. But when I did the transformation with MATLAB's fft I found the result to be of complex numbers because I got nonzero values of the imaginary part of the spectrum. And the magnitude of these imaginary parts is not of the order where you can safely neglect it with respect to the real part. More clearly I will let you have a look at what I have done in respect to your suggestion.

[Matthew Crema]

All,

I've been playing with this a little bit, and I have written some code that works. This is an edited version of code posted earlier, let me know if you want the previous version which was quite flawed. Some comments:

1) This code uses the same number of points in both the time and frequency domain, but I agree with your earlier suggestion that zero-padding could be a good idea. It will make your spectrum "look better" in the frequency domain, but note that as long as Nyquist's criterion is reasonably satisfied it is not necessary to do any zero-padding to reconstruct your time domain waveform. And if there is a possibility that the zero-padding is adding to the confusion, I suggest not doing it.

2) The code uses a different analytic expression that I found on the web.

Chirped Pulse Transform

I did not check the math, but the expression toward the end (on pg 3: "And so finally the electric field spectrum is...") seems to be correct based on the numerical outcome of this code.

3) I ask you to check again that the quantity in one of your exponents is not dimensionless. See comments in code below. Ask yourself what are the units of w0 and w1 (you can figure this out looking at your time domain expression b/c the argument to cosine must also be dimensionless), then check to see what happens if you plug these units into your expression for the Fourier Transform. A side note, it helps (us and you) if you write your expressions in terms of meaningful variables whose units can be checked instead of carrying around the numbers -15 and 7 (and the number 200 which seems to come out of nowhere). And if there are "hidden constants" ( = 1? ) it would help if you included them and told us the units.

3) I got rid of the "time shift and shift back" strategy that I suggested earlier as it just made things more complicated.

Hope this helps. Fun problem anyway.

%%Constants
Npts = 4096;    % Number of points in discrete signals [samples]
Fs = 128;       % Sampling Rate [samples/sec]
Dur = Npts/Fs;  % Duration of time domain signal [sec]
n = (0:Npts-1); % Indexing vector for both time and freq vectors [samples]
time = n/Fs;  % Time vector [sec]
freq = n/Dur;   % Frequency vector [Hz]

% Shift time vector to represent non-causal function and frequency vector
% to represent negative frequencies
halfpt = floor(Npts/2)+1;
% First halves of TSHFT and FSHFT represents negative time/freq, second
% half represents positive time/freq. Use fftshift to put negative halves
% on the right before taking fft or ifft.
tshft = time - time(halfpt);
fshft = freq - freq(halfpt);

%%Compute Chirp-Pulse For Given Set of Parameters

% 1) Define "Chirp" frequency sweep
f0 = -15/(2*pi);    % Instanteneous Freq at beginning of chirp (t=0) [Hz]
w0 = 2*pi*f0;       % Starting Radian Frequency at t = 0 [rad/sec]
k = 7/(2*pi);       % Chirp rate (f1-f0)/Dur [1/sec^2]
w1 = 2*pi*k;        % Chirp rate [rad/sec^2]
% Compute f1 from w0, w1 and Dur
f1 = f0 + k*Dur;  % Instanteneous Freq at end of chirp (t=end) [Hz]

% Compute Chirp frequency and phase as functions of time
f_t = f0 + k*tshft;                 % Instantaneous Frequency [Hz]
Phi_t = 2*pi*(f0 + k*tshft).*tshft;  % Instantaneous Phase [rad]

% 2) Define "Pulse" envelope
Tau = 0.5;  % RMS Width of Gaussian Envelope [sec]
A = 1/(2*Tau)^2;  % Envelope parameter and Real part of z0 [1/sec^2]

% Compute Amplitude of Gaussian Envelope as a function of time
env_t = exp(-A*tshft.^2);  % Envelope [physical units e.g. Voltage]

% 3) Compute Chirped Pulse (these expressions are equivalent expressions,
% choose your favorite). Check that all args to cos and exp are
% dimensionless.
Et_orig = env_t .* cos(Phi_t);
Et_orig1 = env_t .* cos(w0*tshft + w1*tshft.^2);
Et_orig2 = real(exp(-A*tshft.^2-1j*Phi_t));

% 4) Compute fft of Et_orig to later compare with Ew_analytical
Ew_correct = fft(fftshift(Et_orig));

% Real-valued time domain signal obtained by ifft
Et_correct = fftshift(ifft(Ew_correct));

plot(time, Et_orig, time, Et_correct, '--')

%%Build Spectrum from Analytical Exprssion and Compute IFFT
w = 2*pi*fftshift(fshft);  % Radian Frequency Vector [rad/sec]

% Formulation 1 (seems correct)
z0 = A + 1j*w1;  % Substitute into a single complex number [rad/sec^2]
% Check units! All args to exp are dimensionless. Good!
E1 = sqrt(1./(2*conj(z0))) .* exp(-(w-w0).^2./(4*conj(z0))); % Pos Freqs
E2 = sqrt(1./(2*z0)) .* exp(-(w+w0).^2./(4*z0));  % Neg Freqs

% Combine and scale, not clear why scale factor is correct
Ew_analytical = sqrt(2*pi)*Fs/2 * (E1 + E2);

% Alternative formulation (likely incorrect)
N = 200;    % Don't know what this is, assuming it has dimension rad/sec^2
% Note that it is not possible for BOTH of the following terms in the
% exponent to be dimensionless.
% 1) -(w+w0).^2/N. If N has units of rad/sec^2 this IS dimensionless
% 2) -w1*1i*(w+w0).^2/N. Has dimension! [rad/sec^2] No Good!
% A "hidden constants" claim needs to be clarified.
E1_test = exp(-(w+w0).^2/N).*exp(-w1*1i*(w+w0).^2/N);
E2_test = exp(-(w-w0).^2/N).*exp(w1*1i*(w-w0).^2/N);

% Questionable Spectrum of a chirped pulse. Remove "_test" from the name of
% this variable to see how your expression fares.
Ew_analytical_test = -sqrt(2*pi)*Fs/2 * (E1_test + E2_test);

% Real-valued time domain signal obtained by ifft
Et_fromifft = fftshift(ifft(Ew_analytical));

% Discard roundoff error
Et_fromifft = Et_fromifft - ...
    real(Et_fromifft).*(abs(real(Et_fromifft))<1e-3) - ...
    1j*imag(Et_fromifft).*(abs(imag(Et_fromifft))<1e-3);

%%Plots

% Plot Time Waveforms
tax = subplot(3,1,1);
plot(tax, tshft, Et_orig, tshft, Et_fromifft, 'r--', tshft, ...
    Et_orig-Et_fromifft, 'k:')
lg = legend(tax, 'Original', 'Reconstructed', 'Difference');
set(lg, 'FontSize', 6)
ylabel(tax, 'x(t)')
xlabel(tax, 'Time (sec)')
% Determine and set appropriate time axis limits
set(tax, 'Xlim', [-1 1] * min(6*Tau, max(abs(tshft))))
% Plot Real and Imaginary Parts of Fourier Transforms
fax(1) = subplot(3,1,2);
plot(fax(1), fshft, fftshift(real(Ew_correct)), fshft, ...
    fftshift(real(Ew_analytical)), 'r--', fshft, ...
    fftshift(real(Ew_correct-Ew_analytical)), 'k:')
xlabel(fax(1), 'Frequency (Hz)')
ylabel(fax(1), 'Real(FFT(x))')
fax(2) = subplot(3,1,3);
plot(fax(2), fshft, fftshift(imag(Ew_correct)), fshft, ...
    fftshift(imag(Ew_analytical)), 'r--', fshft, ...
    fftshift(imag(Ew_correct-Ew_analytical)), 'k:')
xlabel(fax(2), 'Frequency (Hz)')
ylabel(fax(2), 'Imag(FFT(x))')
% Determine and setappropriate frequeny axis limits
Ew_dB = (20*log10(abs(fftshift(Ew_correct))));
fmax = fshft(find(Ew_dB-max(Ew_dB)>-60, 1, 'last'));
set(fax, 'Xlim', [-1 1]* min(fmax, Fs/2))
linkaxes(fax, 'x')

[Matthew Crema]

And, you can modify my code to test your simpler example by changing three lines to

f0 = -15/(2*pi); becomes f0 = 4/(2*pi);

k = 7/(2*pi); becomes k = 0;

Tau = 0.5, becomes Tau = 2;

Which is another good reason to use variables in your code.

Best,

Matt

Answer 1

[Imam]

I appreciate you efforts Matthew, I tried your last code with the simpler example (no chirp) and it turns out that it gives similar result to the last picture. However if one computes the spectrum in analytical way one should get a real function of the spectrum, i.e. it must not have imaginary part.

And by the way I just found where the problem was. I think for some reason we shouldn't deliberately choose the range of the function to be transformed using either fft or ifft. Precisely, one should not start with a function centered at zero at horizontal axis (in my case the time axis). Instead the starting function should start from zero all the way to the positive region of horizontal axis, only then the application of fft or ifft will yield a correct result. I found this after looking up to the code made by my professor designed for transforming a function that is centered at zero. In his algorithm he first aligned the input function in such a way that its right half was moved to the left, and the left half to the right. After that apply fft and so on, the result matches the manual calculation.

Anyway thanks for the link about chirped pulse transformation. The expression on that paper is the same as what I calculated before, so that now I don't really need to worry if my calculation hold any mistake.

[Matthew Crema]

That's a very good point: the imaginary part of the Fourier Transform for your simple example (which is an even function) should be 0. To use the DFT on non-causal functions like these one must keep in mind that the DFT is periodic and put the "negative" times to the right hand side of the MATLAB array (using fftshift).

I've revised the code in my previous comment, and this should get you closer.

-Matt