matlab使用教程(40)—二维傅里叶变换和多项式插值

1使用 FFT 进行多项式插值

        使用快速傅里叶变换 (FFT) 来估算用于对一组数据进行插值的三角函数多项式的系数。

1.1数学中的 FFT

        FFT 算法通常与信号处理应用相关，但也可以在数学领域更广泛地用作快速计算工具。例如，通常通过解算简单的线性系统来计算用于对一组数据进行插值的 n 次多项式的系数c i：

        在 19 世纪初研究小行星轨道时，卡尔·弗里德里希·高斯发现了一种数学捷径，即通过将问题分解为多个较小的子问题，然后将结果合并来计算多项式插值的系数。他的方法等同于估算其数据的离散傅里叶变换。

1.2 小行星数据插值

        在高斯的论文中，他描述了一种估算小行星智神星轨道的方法。他从以下 12 个二维数据点 x 和 y 开始。

x = 0:30:330;
y = [408 89 -66 10 338 807 1238 1511 1583 1462 1183 804];
plot(x,y,'ro')
xlim([0 360])

        高斯使用以下形式的三角函数多项式为小行星的轨道建模。

        使用 fft 计算多项式的系数。

m = length(y);
n = floor((m+1)/2);
z = fft(y)/m;
a0 = z(1); 
an = 2*real(z(2:n));
a6 = z(n+1);
bn = -2*imag(z(2:n));

        在原始数据点上绘制插值多项式。

hold on
px = 0:0.01:360;
k = 1:length(an);
py = a0 + an*cos(2*pi*k'*px/360) ...
 + bn*sin(2*pi*k'*px/360) ...
 + a6*cos(2*pi*6*px/360); 
plot(px,py)

2二维傅里叶变换

        fft2 函数将二维数据变换为频率空间。例如，您可以变换二维光学掩膜以揭示其衍射模式。

2.1二维傅里叶变换原理

        以下公式定义 m×n 矩阵 X 的离散傅里叶变换 Y。

        计算 X 的二维傅里叶变换等同于首先计算 X 每列的一维变换，然后获取每行结果的一维变换。换言之，命令 fft2(X) 等同于 Y = fft(fft(X).').' 。

2.2二维衍射模式

        在光学领域，傅里叶变换可用于描述平面波入射到带有小孔的光学掩膜上所产生的衍射模式 [1]。本示例对光学掩膜使用 fft2 函数来计算其衍射模式。

        创建用于定义带有小圆孔的光学掩膜的逻辑数组。

n = 2^10; % size of mask
M = zeros(n);
I = 1:n; 
x = I-n/2; % mask x-coordinates 
y = n/2-I; % mask y-coordinates
[X,Y] = meshgrid(x,y); % create 2-D mask grid
R = 10; % aperture radius
A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R
M(A) = 1; % set mask elements inside aperture to 1
imagesc(M) % plot mask
axis image

        使用 fft2 计算掩膜的二维傅里叶变换，并使用 fftshift 函数重新排列输出，从而使零频率分量位于中央。绘制生成的衍射模式频率。蓝色指示较小的幅值，黄色指示较大的幅值。

DP = fftshift(fft2(M));
imagesc(abs(DP))
axis image

        为增强小幅值区域的细节，需绘制衍射模式的二维对数。极小的幅值会受数值舍入误差影响，而矩形网格则会导致径向非对称性。

imagesc(abs(log2(DP)))
axis image