使用“全相位FFT”方法，计算啁啾高斯脉冲的相位谱

前几天，一位好友说用FFT写了一段程序，来计算啁啾高斯脉冲的相位谱，总是得不到正确的结果。后来查找资料发现，用普通FFT算法求相位谱，可能会遇到所得相位谱与实际相位谱差别很大的情况。可以采用“全相位FFT”方法来计算相位谱，这种方法的计算结果比较接近信号的实际相位谱。
啁啾高斯脉冲的时域和频率表达式如图1所示。

               图1 啁啾高斯脉冲的时域和频域表达式，C为啁啾量

使用“全相位FFT”方法求啁啾高斯脉冲相位谱的代码如下(MATLAB)：

%==========================================================================
%Name: apFFT_chirp_Gauss.m，
%Desc: 计算啁啾高斯脉冲的幅度谱和相位谱,使用全相位FFT算法
%Parameter:
%Return:
%Author: yoyoba([e-mail])
%Date: 2013-06-29
%Modify: 2013-06-29
%==========================================================================
close all;
clc;
clear all;
T0=1e-12; %高斯脉冲时间宽度1ps
N=2^11;              %采样点2048
Ts=100*T0/N; %时域取样间隔，时间窗口宽度为100*T0，
n=-N+1:2*N-1;
t=n*Ts; %取样时刻
C=2; %简单一些，直接给出初始啁啾量
g=exp(-0.5*(1+1i*C)*(t/T0).^2); %带有啁啾的高斯脉冲序列，原始脉冲强度默认为1
k=floor(-(N-1)/2:(N-1)/2);
D=1/(N*Ts); %频率分辨率
f=k*D; %频率坐标
%以下代码作用：为了使用全相位FFT算法，需对高斯脉冲序列进行预处理
%=============全相位FFT预处理开始================
win=hanning(N)';
win1=hanning(N)';
win2=conv(win,win1);
win2=win2/sum(win2);
w=pi*2;
g1=g(1:2*N-1);
y1=g1.*win2;
y1a=y1(N:end) + [0 y1(1:N-1)];
%=============全相位FFT预处理结束================
apFFT_g=fftshift(fft(y1a,N)); %FFT变换
figure;
plot(t*1e+12,abs(g).^2/max(abs(g).^2),'linewidth',2);
xlabel('time(ps)');
ylabel('Intensity(a.u.)');
axis([-3,3,0,1]);
title('Time-Intensity');
figure;
plot(t*1e+12,phase(g),'linewidth',2);
xlabel('time(ps)');
ylabel('phase(rad)');
axis([-100,100,-200,700]);
title('Time-Phase');
figure;
plot(f*1e-12,abs(apFFT_g).^2/max(abs(apFFT_g).^2),'linewidth',2);
xlabel('frequency(THz)');
ylabel('Intensity(a.u.)');
axis([-2,2,0,1]);
title('Frequency-Intensity');
figure;
plot(f*1e-12,phase(apFFT_g),'linewidth',2);
xlabel('frequency(THz)');
ylabel('phase(rad)');
axis([-5,5,-70,10]);
title('Frequency-Phase');

运行结果如下：图2~5分别表示时域强度谱、时域相位谱、频域强度谱、频域相位谱。

                  图2 时域强度谱

                图3 时域相位谱

               图4 频域强度谱

               图5 频域相位谱

apFFT_chirp_Gauss.rar