Transformation from spectrum into pulse
21 views (last 30 days)
Show older comments
I have a chirped pulse of the form
E(t) = exp(-t^2)*cos(-15t+7t^2).
Which looks like this
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/145215/image.jpeg)
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,
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/145216/image.jpeg)
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.
8 Comments
Matthew Crema
on 29 Aug 2014
Edited: Matthew Crema
on 30 Aug 2014
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.
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
on 29 Aug 2014
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
Answers (1)
Imam
on 29 Aug 2014
1 Comment
Matthew Crema
on 30 Aug 2014
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
See Also
Categories
Find more on Vibration Analysis in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!