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

前两天，一位好友让我看一段程序，使用FFT算法计算啁啾高斯脉冲的相位谱，总是得不到正确结果。后来查找资料发现，普通FFT算法求相位谱可能会与实际相位谱差别较大。如果改用“全相位FFT”算法，可以得到较为理想的相位谱。啁啾高斯脉冲的时域及频域表达式如图1所示。

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

以下程序使用“全相位FFT算法”，计算其相位谱：

%==========================================================================
%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