此内容理论部分出自书《Numerical Simulation Of Optical Wave Propagation With Example In Matlab》第二章 "Digital Fourier Transforms"

---

1、基本原理：

根据傅里叶变换的定义公式：

         （1）

当对上述公式进行关于x的n阶求导数时，则有：

              （2）

因此，对公式两端取逆傅里叶变换则可得对g(x)公式的对x的n阶求导数的结果。 

---

2、Matlab Code:

• 首先是主函数：

% example_derivative_ft.m
close all;clear all;clc;

N = 64;                 % number of samples
L = 6;                  % grid size [m]
delta = L/N;            % grid spacing [m]
x = (-N/2 : N/2-1) * delta;
w = 3;                  % size of window (or region of interest) [m]
window = rect(x/w);     % window function [m]
g = x.^5 .* window;     % function
% computed derivatives
gp_samp = real(derivative_ft(g, delta, 1)) .* window;
gpp_samp = real(derivative_ft(g, delta, 2)) .* window;
% analytic derivatives
gp = 5 * x.^4 .* window;
gpp = 20 * x.^3 .*  window;

figure,
plot(x,g,x,gp,'-.',x,gp_samp,'x',...
    x,gpp,'--',x,gpp_samp,'+','linewidth',1.2);
xlim([-1,1]); ylim([-20, 20]);
xlabel('x [m]')
legend({'g(x)',"g'(x) analytic","g'(x) numerical", ...
    "g''(x) analytic","g''(x) numerical"},...
    'Location','southeast','FontSize',12);
grid on

Note: 由于example中使用的g(x)函数和他的前几阶到束均不是带限函数，故使用一个window函数（rect(x/w)）来限制显示信号的范围，并减轻计算频谱的混叠。

• 调用子函数实现公式（2）：

function der = derivative_ft(g, delta, n)
 % function der = derivative_ft(g, delta, n)
 % performing a one-dimensional discrete derivative
 N = size(g,2); % number of samples in g
%  grid spacing in the frequency domain
 df = 1/(N*delta);
 fx = (-N/2:N/2-1)*df; % frequency values
 k = (1i*2*pi*fx).^n;
 G = ft(g,delta);
 der = ift(k.*G,df);

逆傅里叶变换： 

function g = ift(G, delta_f)
% function g = ift(G, delta_f)
    g = ifftshift(ifft(ifftshift(G)))...
        * length(G) * delta_f;

傅里叶变换： 

function G = ft(g, delta)
% function G = ft(g, delta)
G = fftshift(fft(fftshift(g)))*delta;

Window function:

function y = rect(x,D)
% matlab code for evaluating the rect function
    if nargin == 1
        D = 1;
    end
    x = abs(x);
    y = double(x<D/2);
    y(x == D/2) = 0.5;

---

3、运行结果：