Chapter 2 - 利用DFT求函数导数
作者:互联网
此内容理论部分出自书《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、运行结果:
标签:Chapter,function,求函数,ft,gp,DFT,window,delta,derivative 来源: https://blog.csdn.net/qq_33801693/article/details/118932815