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

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?

Answer 1

[Wayne King]

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.

[Dong-Gyu Jang]

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.

[Wayne King]

It is not just MATLAB's problem.

Answer 2

[Dong-Gyu Jang]

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.

[TheMellow]

This is very helpful, thanks! Not sure whether you meant to do this but there is a bracketing error in your code. This is the fixed version:

y(1:mid_t_s) = exp(-1.38629*(t_s(1:mid_t_s)/t_r).^2).*cos(w*t_s(1:mid_t_s));

y(mid_t_s:length(t_s)) = exp(-1.38629*(t_s(mid_t_s:end)/t_f).^2).*cos(w*t_s(mid_t_s:end));

[Jonathan]

The Fig. 1 plotted using this code looks strange.

[Dong-Gyu Jang]

Thanks.

I attached the code right now and believe the code in my answer fixed the problem what you mentioned.