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

15 次查看（过去 30 天）
Dong-Gyu Jang
Dong-Gyu Jang 2012-8-5
评论： Dong-Gyu Jang 2016-2-3
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?

请先登录，再回答此问题。

采纳的回答

Wayne King
Wayne King 2012-8-5
编辑：Wayne King 2012-8-5
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 个评论
Dong-Gyu Jang
Dong-Gyu Jang 2012-8-5
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.

请先登录，再进行评论。

更多回答（1 个）

Dong-Gyu Jang
Dong-Gyu Jang 2012-8-16
编辑：Dong-Gyu Jang 2016-2-3
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 个评论
显示 1更早的评论
Dong-Gyu Jang
Dong-Gyu Jang 2016-2-3
Thanks.
I attached the code right now and believe the code in my answer fixed the problem what you mentioned.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Estimation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by