Zabbix告警消息推送至kafka

JSON技术

  返回  

Chapter 2 - 利用DFT求函数梯度

2021/7/20 14:09:33 浏览:

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


1、基本原理:

现结利用DFT求函数导数的内容,Chapter 2 - 利用DFT求函数导数_阿甘番外篇的博客-CSDN博客,利用DFT求函数的梯度。

针对二维函数,分别对函数的两个自变量求n阶偏导数,可得:

\Gamma \left \{ \frac{\partial ^{n}}{\partial x^n} g(x,y)\right \} = (i2\pi f_x)^n\Gamma\left \{ g(x,y) \right \}          (1)

\Gamma \left \{ \frac{\partial ^{n}}{\partial y^n} g(x,y)\right \} = (i2\pi f_y)^n\Gamma\left \{ g(x,y) \right \}         (2)

 则此函数的梯度在n=1的情况下表示为:

\bigtriangledown g(x,y) = \Gamma ^{-1}\left \{ i2\pi f_x\Gamma \left \{ g(x,y) \right \} \right \}\hat{i}+ \Gamma ^{-1}\left \{ i2\pi f_y\Gamma \left \{ g(x,y) \right \} \right \}\hat{j}       (3)

因此,只需要利用DFT求解函数的(逆)傅里叶变换,即可获得函数g(x)的梯度 


  2、Matlab Code:

  • 主函数:
% example_gradient_ft.m
close all;clear all;clc;

N = 64; % number of samples
L = 6; % grid size [m]
delta = L/N; % grid spacing [m]
xx = (-N/2:N/2-1)*delta;
[x,y] = meshgrid(xx);
g = exp(-(x.^2+y.^2)); % signal function
% compute derivatives
[gx_samp, gy_samp] = gradient_ft(g, delta);
gx_samp = real(gx_samp);
gy_samp = real(gy_samp);
% analytic derivatives
gx = -2*x.*exp(-(x.^2+y.^2));
gy = -2*y.*exp(-(x.^2+y.^2));

figure,
imagesc(xx,xx,g); colorbar;
xlabel('x [m]'); ylabel('y [m]');
title('Original Function g(x,y)')
axis([-1 1 -1 1]);axis square

figure
quiver(x, y, gx_samp, gy_samp)
xlabel('x [m]'); ylabel('y [m]');
title('Numeric \nabla g(x,y)')
axis([-1 1 -1 1]);axis square

figure
quiver(x, y, gx, gy)
xlabel('x [m]'); ylabel('y [m]');
title('Analytic \nabla g(x,y)')
axis([-1 1 -1 1]);axis square
  • 调用子函数:
function [gx, gy] = gradient_ft(g, delta)
% function [gx gy] = gradient_ft(g, delta)
% computing the discrete gradient of a function using FTs.

    N = size(g,1); % number of samples per side in g
    % grid spacing in the frequency domain
    df = 1/(N*delta);
    fx = (-N/2 : N/2-1)*df;
    [fx, fy] = meshgrid(fx);
    G = ft2(g, delta);
    gx = ift2(1j*2*pi*fx .* G, df);
    gy = ift2(1j*2*pi*fy .* G, df);

二维傅里叶变换函数:

function G = ft2(g, delta)
% function G = ft2(g, delta)
    G = fftshift(fft2(fftshift(g))) * delta^2;

二维逆傅里叶变换函数:

function g = ift2(G, delta_f)
% function g = ift2(G, delta_f)
    N = size(G,1);
    g = ifftshift(ifft2(ifftshift(G))) * (N * delta_f)^2;

运行结果:

 

联系我们

如果您对我们的服务有兴趣,请及时和我们联系!

服务热线:18288888888
座机:18288888888
传真:
邮箱:888888@qq.com
地址:郑州市文化路红专路93号