Why does spectral phase from the fft of gaussian pulse shows sawteeth?

27 views (last 30 days)
Hello.
I'm interested in spectral analysis on laser pulse. I made a test code in which gaussian pulse is analyzed in FFT. The analytic answer is that the spectral amplitude is also gaussian profile while spectral phase is all zero.
clear; close all;
% Temporal profile
delta_t = 0.1; % Sampling interval.
t = -90:delta_t:90;
t_p = 35; wavelength = 790E-9; c = 3E8; k = 2*pi/wavelength; w = k*c; w = w*1E-15; % It is Just for getting normalized central frequency. Not imporant for this code.
E_field = exp(-1.38629*(t/t_p).^2).*cos(w*t);
figure(1); subplot(2,2,1); plot(t,E_field); title('Temporal envelope'); xlabel('Time(fs)'); ylabel('E-field');
% Spectral profile
N = 10*length(t); % number of FFT data points is 10 times data points of temporal signal.
spectrum = fft(E_field,N);
frequency = (0:N-1)/(N*delta_t);
figure(1); subplot(2,2,2); plot(frequency,abs(spectrum)); title('Spectrum'); xlabel('frequency(1*10^1^5 Hz)');
% Detailed spectral profile including phase
[pkv,pkl] = max(abs(spectrum(1:fix(N/2)))); % pkl is location of global maxima.
data_width = 50; % Data width over which a left peak is showed in details.
spectral_amplitude = abs(spectrum);
spectral_phase = angle(spectrum);
figure(1); subplot(2,2,3); plot(frequency(pkl - data_width:pkl + data_width),spectral_amplitude(pkl - data_width:pkl + data_width)); title('Detailed spectrum'); xlabel('frequency(1*10^1^5 Hz)'); figure(1);
subplot(2,2,4); plot(frequency(pkl - data_width:pkl + data_width),spectral_phase(pkl - data_width:pkl + data_width)); title('Spectral phase'); xlabel('frequency(1*10^1^5 Hz)'); phase = spectral_phase(pkl)
Number of FFT data points N are far larger than that of temporal signal and is also integral multiple of delta_t (= 1/sampling frequency) to avoid so called leakage problem.
It is found that the amplitude is really good however, phase profile is sawteeth shape far from what I expected.
I believe that there is no reason why phase continuously decreases thus phase jumps are occured.
Is it due to some numerical error? what kind of error should I consider? How can I solve this problem thus exact phase analysis is being made?

Accepted Answer

Wayne King
Wayne King on 5 Aug 2012
Edited: Wayne King on 5 Aug 2012
I don't think there is a problem, I think if you plot the following
subplot(211)
plot(frequency,abs(spectrum))
subplot(212)
plot(frequency,unwrap(angle(spectrum)))
You see that the phase is relatively constant except in the location of the positive and negative frequency peaks of the sinusoidal component.
You have to keep in mind that most of your DFT coefficients have small magnitudes and small variations in those values can cause jumps in the phase. Try a simple experiment with a single unwindowed cosine
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t);
xdft = fft(x);
plot(abs(xdft))
figure
plot(angle(xdft))
The phase is making apparent jumps throughout, but you have to consider whether those apparent phase jumps really matter. The truth is that most of them are associated with DFT coefficients with magnitudes on the order of 10^{-13} or so.
This is just the reality of phase estimation with computer implementations of the DFT.
  2 Comments
Dong-Gyu Jang
Dong-Gyu Jang on 5 Aug 2012
Thank you for answering my question as quickly as possible.
Investigation of spectral phase of various laser pulse is really important for making pulse shaping which is my interests. Thus I really need to get spectral phase as exact as possible.
Is it just Matlab's problem or intrictic propert DFT should have?
Is there any other way to get such a good phase? Can I make some reference and subtract that from data thus good phase can be obtainable?
If you know any method or document about this issue, Please let me know that.

Sign in to comment.

More Answers (1)

Dong-Gyu Jang
Dong-Gyu Jang on 16 Aug 2012
Edited: Dong-Gyu Jang on 3 Feb 2016
I finally got right answer!
Physically, gaussian pulse is located around t = 0 thus, symmetry around origin is fullfiled.
However, in DFT view, sample data is not symmetric - data is not N-points circularly even.
Data which is N-points circularly even satisfies x(n) = x(N-n), n = 1,2,3,...N where x(n) n-th sample data.
When N-points circularly even is satisfied and data is real (not complex value), than spectral phase becomes zero.
One who want to know property of DFT in details might need to read http://web.eecs.umich.edu/~fessler/course/451/l/pdf/c5.pdf or http://dsp-book.narod.ru/TDCH/CH-02.PDF
Anyway, following code is modified version and analytically expected zero phase over spectrum is obtained.
The most important modification is to add t = ifftshift(t); before making sample data. By adding this code, sample data satisfies N-points circularly even.
I attached modified M-scripts showing right result.
  3 Comments
Dong-Gyu Jang
Dong-Gyu Jang on 3 Feb 2016
Thanks.
I attached the code right now and believe the code in my answer fixed the problem what you mentioned.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!