首页 > 其他分享> 文章详细

Chapter 2 - 利用DFT求函数导数

2021-07-20 14:02:31 阅读：256 来源： 互联网

标签：Chapter function 求函数 ft gp DFT window delta derivative

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

1、基本原理：

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

$G(f_{x})=\Gamma {g(x)} = \int_{\infty }^{\infty}g(x)e^{-i2\pi f_{x}x}dx$ （1）

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

$\Gamma \left \{ {\frac{d^{n}}{dx^{n}}}g(x) \right \}=(i2\pi f_x)^n\Gamma\left \{ g(x) \right \}$ （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、运行结果：

标签：Chapter,function,求函数,ft,gp,DFT,window,delta,derivative
来源： https://blog.csdn.net/qq_33801693/article/details/118932815

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

Chapter 2 - 利用DFT求函数导数