标签: dft

我试图通过一个例子验证我对Numpy的FFT的理解:傅立叶变换exp(-pi*t^2)应该是exp(-pi*f^2)在直接变换上没有应用缩放时.

但是,我发现为了获得这个结果,我需要将FFT的结果乘以一个因子dt,即我函数上两个采样点之间的时间间隔.我不明白为什么.有人可以帮忙吗？

这是一个示例代码:

# create data
N = 4097
T = 100.0
t = linspace(-T/2,T/2,N)
f = exp(-pi*t**2)

# perform FT and multiply by dt
dt = t[1]-t[0]
ft = fft(f)  * dt      
freq = fftfreq( N, dt )
freq = freq[:N/2+1]

# plot results
plot(freq,abs(ft[:N/2+1]),'o')
plot(freq,exp(-pi*freq**2),'r')
legend(('numpy fft * dt', 'exact solution'),loc='upper right')
xlabel('f')
ylabel('amplitude')
xlim(0,1.4)

我正在尝试解决python + numpy的问题,其中我有一些类型的函数 我需要与另一个函数进行卷积 .为了优化代码,我执行了f和g的fft,我将它们相乘,然后我执行逆变换以获得结果.

作为进一步的优化,我意识到,由于移位定理,我可以简单地计算一次f(x,y,z)的fft,然后将其乘以相位因子,这取决于 获得fft .特别是,,其中N是x和y的长度.

我尝试用python + numpy来实现这个简单的公式,但由于某些原因我目前无法理解它失败了,所以我要求SO社区的帮助来弄清楚我缺少的东西.

我还提供了一个简单的例子.

In [1]: import numpy as np
In [2]: x = np.arange(-10, 11)
In [3]: base = np.fft.fft(np.cos(x))
In [4]: shifted = np.fft.fft(np.cos(x-1))
In [5]: w = np.fft.fftfreq(x.size)
In [6]: phase = np.exp(-2*np.pi*1.0j*w/x.size)
In [7]: test = phase * base
In [8]: (test == shifted).all()
Out[8]: False
In [9]: shifted/base
Out[9]:
array([ 0.54030231 -0.j        ,  0.54030231 -0.23216322j,
        0.54030231 -0.47512034j,  0.54030231 -0.7417705j ,
        0.54030231 -1.05016033j,  0.54030231 -1.42919168j,
        0.54030231 -1.931478j …

假设信号对应于一年(365天)内的日值.它由全零组成,除了一些稀疏值,它们对应于所有以相同间隔(30天)隔开的孤立峰值.我使用快速傅里叶变换功能获得频谱.

1. 如何摆脱高0Hz的峰值？ 编辑:这与信号的非零均值性质有关.有关详细信息,请参阅此帖子.

2. 然后第一个峰值为12Hz,这在某种程度上是预期的.然而,峰值也存在于24Hz,36Hz,48Hz ...... 这是一个走样问题吗？如何摆脱它？

以下是我的代码.它在Octave中进行了测试,但它也应该在Matlab中运行

close all
clear all

T = 1/365; % period
samp_freq = 1/T;  % sample frequency
t=0:T:2; % overall time span is two years

% build signal
x= zeros(length(t),1);    
for i=1:length(t)
    if mod(i,30) == 0
       x(i) = 100; 
    else
        x(i) = 0;
    end
end

figure(1)
plot(t,x)
grid
xlabel("Time [years]")
ylabel("Signal amplitude")


y=fft(x);
N = length(x);

for i=1:N
    f(i) = (i-1)*samp_freq/N;
end

figure(2)
plot(f,abs(y))
xlabel("Frequency")
ylabel("Signal amplitude")

figure(3)
plot(f(1:80),abs(y(1:80)))
xlabel("Frequency")
ylabel("Signal amplitude")

我正在尝试实现一个基于DFT的8波段均衡器,仅用于学习.为了证明我的DFT实现工作,我输入了一个音频信号,对其进行了分析,然后再次重新合成它,而没有对频谱进行任何修改.到现在为止还挺好.

我正在使用所谓的"计算DFT的标准方式",这是通过相关性.该方法计算长度为N/2 + 1个样本的实部和虚部.为了减弱我正在做的频率:

float atnFactor = 0.6;
Re[k] *= atnFactor;
Im[k] *= atnFactor;

其中'k'是0到N/2范围内的指数,但我在再合成后得到的是一个轻微的失真信号,特别是在低频时.

输入信号采样率是44.1 khz,因为我只想要一个8波段均衡器,我一次只能给DFT 16采样,所以我有8个频率区可供使用.

有人能告诉我我做错了什么吗？我试图在互联网上找到关于这个主题的信息,但找不到任何信息.

提前致谢.

我是OpenCV和图像处理算法的新手.我需要在C++中的OpenCV中进行逆离散傅里叶变换,但我不知道如何.我在互联网上搜索,但没有找到答案.我正在使用此页面中的代码在我的程序中进行傅里叶变换:http://opencv.itseez.com/doc/tutorials/core/discrete_fourier_transform/discrete_fourier_transform.html.我试图对该代码做反向,但我不知道我在哪里做错了.我的代码在这里(我认为整个代码都是错误的):

void doFourierInverse(const Mat &src, Mat &dst) {
  normalize(src, dst, 0, -1, CV_MINMAX); 

  int cx = dst.cols/2;
  int cy = dst.rows/2;

  Mat q0(dst, Rect(0, 0, cx, cy));   
  Mat q1(dst, Rect(cx, 0, cx, cy));  
  Mat q2(dst, Rect(0, cy, cx, cy));  
  Mat q3(dst, Rect(cx, cy, cx, cy)); 

  Mat tmp;         
  q0.copyTo(tmp);
  q3.copyTo(q0);
  tmp.copyTo(q3);

  q1.copyTo(tmp);      
  q2.copyTo(q1);
  tmp.copyTo(q2);

  dst = dst(Rect(0, 0, dst.cols & -2, dst.rows & -2));

  exp(dst, dst);
  dst -= Scalar::all(1);  

  Mat planes[2];

  polarToCart(dst, Mat::zeros(dst.rows, dst.cols, dst.type()), planes[0], planes[1]); …

我尝试计算这个数组的DFT x_1.它必须简单,但我的价值太低了.我的代码出了什么问题？

请不要链接到其他示例 - 只是寻找我自己的代码的修复程序.

#include <iostream>
#include <complex>
#include <cassert>

int main ()
{
    const unsigned int N = 20;

    const double x_1[N] = {0, 0.3, 0.6, 0.8, 1, 1, 0.9, 0.7, 0.5, 0.2, 0.2, 0.5, 0.7, 0.9, 1, 1, 0.8, 0.6, 0.3, 0};

    for(unsigned int k = 0; k < N; k++)
    {
        std::complex<double> sum(0.0,0.0);
        for(unsigned int j = 0; j < N; j++)
        {
            int integers = -2*j*k;
            std::complex<double> my_exponent(0.0, M_PI/N*(double)integers);
            sum += x_1[j] …

根据卷积定理，时域中的卷积是fft域中的乘积。通过正确的零填充，它可以工作：

% convolution in time domain
a = [1 2 3];
b = [4 5 6];
c = conv(a,b);

a_padded=[a 0 0]; b_padded=[b 0 0];
c_bis=ifft(fft(a_padded).*fft(b_padded));
% we do find c_bis=c

然而，这个定理也应该以相反的方式工作，时域中的乘积是 fft 域中的卷积。我不明白这部分：

d = a.*b;
D=conv(fft(a_padded),fft(b_padded));
d_bis=ifft(D);

这给出了 d_bis 的复向量。如何使用频域中的卷积来反转时域中的逐点乘积？

我目前正在尝试使用Kiss FFT将FFT实施到AVR32微控制器中，以进行信号处理。而且我的输出有一个奇怪的问题。基本上，我会将ADC样本（使用函数发生器进行测试）传递到fft（真实输入，256 n大小）中，并且检索到的输出对我来说很有意义。但是，如果我将汉明窗应用于ADC样本，然后将其传递给FFT，则峰值幅度的频率仓是错误的（并且与之前没有开窗的结果不同）。ADC样本具有DC偏移，因此我消除了偏移，但仍不适用于窗口样本。

以下是通过rs485的前几个输出值。第一列是没有窗口的fft输出，而第二列是有窗口的输出。从第1列开始，峰值位于第6行（6 x fs（10.5kHz）/ 0.5N）为我提供了正确的输入频率结果，其中第2列在第2行具有峰值幅度（直流仓除外），这对我来说没有意义。任何建议都会有所帮助。提前致谢。

    488260 //直流仓
    5 97
    5 41
    5 29  
    4 26
    10 35
    133 76
    33 28
    21 6
    17 3

kiss_fft_scalar zero;
memset(&zero,0,sizeof(zero));
kiss_fft_cpx fft_input[n];
kiss_fft_cpx fft_output[n];

for(ctr=0; ctr<n; ctr++)
{
    fft_input[ctr].r = zero;
    fft_input[ctr].i = zero;
    fft_output[ctr].r =zero;
    fft_output[ctr].i = zero;
}

// IIR filter calculation

for (ctr=0; ctr<n; ctr++)
{       
    // filter calculation
    y[ctr] = num_coef[0]*x[ctr];

    y[ctr] += (num_coef[1]*x[ctr-1]) - (den_coef[1]*y[ctr-1]);
    y[ctr] += (num_coef[2]*x[ctr-2]) - (den_coef[2]*y[ctr-2]);
    //y1[ctr] += y[ctr] - 500; …

我的目标DFT是在 OpenCV 中获取图像。

使用dft函数，我可以计算它，然后通过计算它的大小来绘制它（然后，应用对数并最终对其进行归一化以绘制 0 到 1 之间的值）。

我的结果是，对于下面的图像，我向您展示的结果（交换以在图像的中心具有较低的频率）：

但是，如果我将它与使用其他工具（如Halcon ）获得的结果进行比较，这对我来说似乎不正确，因为它似乎具有非常“高”的值（我的意思是 OpenCV DFT 量级）：

我认为可能是以下原因：

1. DFT（在OpenCV）和FFT（Halcon）的区别
2. 我正在执行的操作是为了显示OpenCV 中的大小。

第一个有我很难分析的问题，OpenCV没有FFT功能，Halcon也没有DFT功能（当然如果我没有错的话），所以我可以不直接比较。

第二个是我工作时间最长的地方，但我仍然找不到原因。

img这是我用来绘制大小的代码（这是我的 DFT 图像）：

// 1.- To split the image in Re | Im values
Mat planes[] = {Mat_<float>(img), Mat::zeros(img.size(), CV_32F)};

// 2.- To magnitude + phase
split(img, planes);

// Calculate magnitude. I overwrite it, I know, but this is inside a function so it …

我需要在图像上应用简单的低通滤波器,但这需要在频域中完成.然而,最终结果是一个非常嘈杂的图像,表明我的过程在某处是不正确的.

我的思想轨迹一直是

1)使图像居中(Testpatt)并应用2D FFT

tpa1 = size(Testpatt, 1); tpa2 = size(Testpatt, 2);
dTestpatt = double(Testpatt);


for i = 1:tpa1
    for j = 1:tpa2

        dTestpatt(i,j) = (dTestpatt(i,j))*(-1)^(i + j);
    end
end

fftdTestpatt = fft2(dTestpatt);

2)生成低通滤波器并将其与傅里叶变换图像相乘(注意:低通滤波器需要在半径范围内为1的圆圈)

lowpfilter = zeros(tpa1, tpa2);
Do = 120;

for i = 1:tpa1
    for j = 1:tpa2


       if sqrt((i - tpa1/2)^2 + (j - tpa2/2)^2) <= Do
           lowpfilter(i,j) = 1;
       end

   end
end


filteredmultip = lowpfilter*fftdTestpatt;

3)采用逆2D FFT的实部,并恢复对中.

filteredimagecent = real(ifft2(filteredmultip));


for i = 1:tpa1
    for j = …