- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我正在尝试构建二维框函数的 FFT 和二维高斯函数的 FFT 的乘积。之后,我找到了逆 FFT,结果我期待这两个函数的卷积。但是,我得到了一个奇怪的单方面结果,如下所示。结果显示在 subplot
的右下角。
我为重现上述 subplot
而编写的 Octave 代码以及我为构建卷积而执行的计算如下所示。谁能告诉我我做错了什么?
clear all;
clc;
close all;
% domain on each side is 0-9
L = 10;
% num subdivisions
N = 32;
delta=L/N;
sigma = 0.5;
% get the domain ready
[x,y] = meshgrid((0:N-1)*delta);
% since domain ranges from 0-(N-1) on both sdes
% we need to take the average
xAvg = sum(x(1, :))/length(x(1,:));
yAvg = sum(y(:, 1))/length(x(:,1));
% gaussian
gssn = exp(- ((x - xAvg) .^ 2 + (y - yAvg) .^ 2) ./ (2*sigma^2));
function ret = boxImpulse(a,b)
n = 32;
L=10;
delta = L/n;
nL = ((n-1)/2-3)*delta;
nU = ((n-1)/2+3)*delta;
if ((a >= nL) && (a <= nU) && ( b >= nL) && (b <= nU) )
ret=1;
else
ret=0;
end
ret;
endfunction
boxResponse = arrayfun(@boxImpulse, x, y);
subplot(2,2,1);mesh(x,y,gssn); title("gaussian fun");
subplot(2,2,2);mesh(x,y, abs(fft2(gssn)) .^2); title("fft of gaussian");
subplot(2,2,3);mesh(x,y,boxResponse); title("box fun");
inv_of_product_of_ffts = abs(ifft2(fft2(boxResponse) * fft2(gssn))) .^2 ;
subplot(2,2,4);mesh(x,y,inv_of_product_of_ffts); title("inv of product of fft");
最佳答案
让我们首先解决您最明显的错误。这是您在频域中计算卷积的方式:
inv_of_product_of_ffts = abs(ifft2(fft2(boxResponse) * fft2(gssn))) .^2 ;
第一个问题是您使用的是*
,即矩阵乘法。在频域中,element-wise乘法相当于频域中的卷积,所以需要用.*
代替。第二个问题是您也不需要 .^2
项和 abs
操作。您正在执行卷积,然后找到绝对值并对每个项进行平方。这不是必需的。删除 abs
和 .^2
操作。
因此,你需要的是:
inv_of_product_of_ffts = ifft2(fft2(boxResponse).*fft2(gssn)));
但是,您将要得到的是这个结果。让我们把它放在一个新图中而不是 subplot
*:
figure;
mesh(x,y,ifft2(fft2(boxResponse).*fft2(gssn)));
title('Convolution... not right though');
您可以看到这是正确的结果……但它没有居中……这是为什么?这实际上是在频域中计算卷积时最常见的问题之一。事实上,即使是最有经验的人也会遇到这个问题,因为他们不了解 FFT 工作原理的内部结构。
这是 MATLAB 和 Octave 如何计算 2D FFT 的结果。具体来说,MATLAB 和 Octave 定义了一个 meshgrid
坐标,水平坐标从 0,1,...M-1
到 0,1,.. .N-1
表示垂直,给我们一个 M x N
结果。这意味着原点/直流分量位于矩阵的左上角,而不是我们通常定义的中心。具体来说,传统的 2D FFT 定义坐标从 -(M-1)/2, ..., (M-1)/2
用于水平和 -(N-1)/2, ..., (N-1)/2
垂直。
引用上述内容,您将信号定义为中心,假设它是原点而不是左上角。为了弥补这一点,您需要添加一个 fftshift
(MATLAB doc,Octave doc),以便 FFT 的输出现在以原点为中心,而不是左上角。回到原来的样子,因此你真的需要:
figure;
mesh(x,y,fftshift(ifft2(fft2(boxResponse).*fft2(gssn))));
title('Convolution... now it is right');
如果您想仔细检查我们的结果是否正确,您可以在空间域中执行直接卷积,我们可以在频域中比较两种方法的结果。
这就是您直接在空间域中计算结果的方式:
figure;
mesh(x,y,conv2(gssn,boxResponse,'same'));
title('Convolution... spatial domain');
conv2
( MATLAB doc , Octave doc ) 在两个信号之间执行二维卷积,'same'
标志确保输出大小是最大的两个信号,可以是高斯滤波器或箱式滤波器。您会看到它是同一条曲线,为简洁起见,我不会在此处显示。
但是,我们可以比较这两个结果,看看它们在元素方面是否相同。执行此操作的一种方法是确定减去结果中的每个元素是否小于某个阈值……比如说……1e-10
:
>> out1 = conv2(boxResponse, gssn, 'same');
>> out2 = fftshift(ifft2(fft2(boxResponse).*fft2(gssn)));
>> all(abs(out1(:)-out2(:)) < 1e-10)
ans =
1
这意味着这两者确实是一样的。
希望这是有道理的!请记住,由于您是在原点位于中心的位置定义信号,一旦找到逆向 FFT,就必须移动,以便原点现在位于中心,而不是顶部 -左角。
*:小注 - 此答案中生成的所有图都是使用 MATLAB R2015a 生成的。但是,该代码在 Octave 中完全有效 - 使用 Octave 4.0.0 进行了测试。
关于matlab - 通过 FFT 执行 2D 卷积时的奇怪行为,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/34081933/
FFT 库(例如 FFTW 或 numpy.fft)通常提供两个函数 fft() 和 ifft()(及其用于实值输入的特殊版本)。这些功能似乎被定义为 ifft(fft(X)) == X 和 fft(
如果我有一个特定大小 M(2 的幂)的 FFT 实现,我如何计算一组大小 P=k*M 的 FFT,其中 k 也是 2 的幂? #define M 256 #define P 1024 comple
下午好! 我正在尝试基于我已有的简单递归 FFT 实现来开发 NTT 算法。 考虑以下代码(coefficients'的长度,让它为m,是2的精确幂): /// /// Calculates the
我正在分析时间序列数据,并希望提取 5 个主要频率分量并将其用作训练机器学习模型的特征。我的数据集是 921 x 10080 。每行是一个时间序列,总共有 921 个。 在探索可能的方法时,我遇到了各
我找不到任何官方文档来证明 scipy.fft 实际上是 numpy.fft.fftpack.fft 的链接。这是显示链接的 iPython session : In [1]: import scip
文档说 np.fft.fft 这样做: Compute the one-dimensional discrete Fourier Transform. 和 np.fft.rfft 这样做: Compu
近一个月来,我一直在与一个非常奇怪的错误作斗争。问你们是我最后的希望。我用 C 编写了一个程序,它集成了 2d Cahn–Hilliard equation在傅里叶(或倒数)空间中使用隐式欧拉 (IE
我一直在制作一个例程,使用 NumPy/Scipy 测量两个光谱之间的相位差。 我已经有了Matlab写的例程,所以我基本上是用NumPy重新实现了函数和相应的单元测试。但是,我发现单元测试失败了,因
我正在研究使用 Renderscript 对大型复杂输入数组执行 FFT。 FFT 是相当标准的,因为它涉及三个循环,但内部循环执行 FFT 中的蝶形运算。因为每个蝴蝶使用数组的不同部分,所以没有明显
我需要通过修改 FFT 结果来均衡音乐样本。 我知道如何获得每个输出虚数的频率,问题是修改这个值以获得“均衡器效果”。 我需要知道如何缩放这个值。 条目大小为 4096 个样本,采样率为 44100
我将在 kiss-fft 之前制定几个计划同时(平行),我可以这样做吗,或者换句话说,kiss-fft 线程安全吗? 谢谢 最佳答案 自述文件: No static data is used. Th
要在频域中插入信号,可以在时域中填充零并执行 FFT。 假设给定向量 X 中的元素数为 N 并且 Y 与 X 相同但在一侧用 N 零填充。然后下面给出相同的结果。 $$\hat{x}(k)=\sum_
我通过相关了解了 DFT 的工作原理,并将其用作理解 FFT 结果的基础。如果我有一个以 44.1kHz 采样的离散信号,那么这意味着如果我要获取 1 秒的数据,我将有 44,100 个样本。为了对其
有人知道 Mayer FFT 的实现吗(我不必花很多时间研究代码)? 我正在尝试执行卷积,ifft 似乎产生了我称之为“镜像”的输出。换句话说,我的内核+信号长度被限制为 N/2 并且占据 n=0..
有人知道 Mayer FFT 的实现吗(我不必花很多时间研究代码)? 我正在尝试执行卷积,ifft 似乎产生了我称之为“镜像”的输出。换句话说,我的内核+信号长度被限制为 N/2 并且占据 n=0..
我有以下代码...请注意#生成正弦曲线下的两行。一个使用比另一个更高的 2pi 精度值,但它们仍然应该给出几乎相同的结果。 import numpy as np import matplotlib.p
我正在努力确保 FFTW 做我认为它应该做的事情,但我遇到了问题。我正在使用 OpenCV 的 cv::Mat。我制作了一个测试程序,给定一个 Mat f,计算 ifft(fft(f)) 并将结果与
我是从事电信项目的计算机程序员。 在我们的项目中,我必须将一系列复数更改为它们的傅立叶变换。因此我需要一个高效的 FFT 代码来满足 C89 标准。 我正在使用以下代码,它运行良好: shor
我目前正在尝试了解 numpy 的 fft 函数。为此,我测试了以下假设: 我有两个函数,f(x) = x^2 和 g(x) = f'(x) = 2*x。根据傅立叶变换定律和 wolfram alph
我一直在使用 FFT,目前正在尝试使用 FFT 从文件中获取声音波形(最终对其进行修改),然后将修改后的波形输出回文件。我得到了声波的 FFT,然后对其使用了反 FFT 函数,但输出文件听起来一点也不
我是一名优秀的程序员,十分优秀!