标签: fftw

我有一个以列主要(Fortran风格)格式存储的2D数据库,我想对每一行进行FFT .我想避免转置数组(它不是正方形).例如,我的数组

fftw_complex* data = new fftw_complex[21*256];

包含条目[r0_val0, r1_val0,..., r20_val0, r0_val1,...,r20_val255].

我可以用它fftw_plan_many_dft来制定一个计划来解决data数组中的21个FFT中的每个FFT,如果它是行 -主要的,例如[r0_val0, r0_val1,..., r0_val255, r1_val0,...,r20_val255]:

int main() {
  int N = 256;
  int howmany = 21;
  fftw_complex* data = new fftw_complex[N*howmany];
  fftw_plan p;

  // this plan is OK
  p = fftw_plan_many_dft(1,&N,howmany,data,NULL,1,N,data,NULL,1,N,FFTW_FORWARD,FFTW_MEASURE);
  // do stuff...

  return 0;
}

根据文档(FFTW手册第4.4.1节),该功能的签名是

fftw_plan fftw_plan_many_dft(int rank, const int *n, int howmany,
                              fftw_complex *in, const int *inembed,
                              int istride, int idist, …

我正在尝试开发一个简单的C应用程序,它可以在WAV文件中的给定时间戳下在特定频率范围内给出0-100的值.

示例:我的频率范围为44.1kHz(典型的MP3文件),我想将该范围分成n个范围(从0开始).然后我需要得到每个范围的幅度,从0到100.

到目前为止我管理的内容:

使用libsndfile我现在能够读取WAV文件的数据.

infile = sf_open(argv [1], SFM_READ, &sfinfo);

float samples[sfinfo.frames];

sf_read_float(infile, samples, 1);

但是,我对FFT的理解相当有限.但我知道为了使振幅达到我需要的范围是必需的.但是我该如何继续前进呢？我找到了FFTW-3库,它似乎适用于此目的.

我在这里找到了一些帮助:https://stackoverflow.com/a/4371627/1141483

并在这里查看了FFTW教程:http://www.fftw.org/fftw2_doc/fftw_2.html

但由于我不确定FFTW的行为,我不知道从这里开始.

另一个问题,假设您使用libsndfile:如果强制读取单引导(使用立体声文件)然后读取样本.那么你真的只会阅读总文件的一半样本吗？其中一半来自频道1,还是自动过滤出来？

非常感谢您的帮助.

编辑:我的代码可以在这里看到:

double blackman_harris(int n, int N){
double a0, a1, a2, a3, seg1, seg2, seg3, w_n;
a0 = 0.35875;
a1 = 0.48829;
a2 = 0.14128;
a3 = 0.01168;

seg1 = a1 * (double) cos( ((double) 2 * (double) M_PI * (double) n) / ((double) N - (double) 1) );
seg2 = a2 * (double) cos( …

我有一个以行主顺序存储的矩阵.我正在尝试使用FFTW计算ub矩阵的DCT,我得到了无意义.在下面的段落中,我将描述问题和我的解决方案,我希望你帮助理解为什么它不起作用.

鉴于一些i和l,我想计算一个子矩阵包括所有行的DCT k为此k mod l == i.例如,假设l = 3和i = 2.在下面的矩阵中,我希望变换的子矩阵用红色标记(2 mod 3 = 2,5 mod 3 = 2,8 mod 3 = 2).

源和目标数组具有相同的布局,转换后的矩阵应存储在目标数组中的相同位置.

void transform(double* src, double* dest, size_t rows, size_t cols, size_t l, size_t i)
{
    int rank = 2;
    fftw_iodim64 dims[] = {
        { rows / l, l, l },
        { cols, rows, rows } };
    fftw_r2r_kind kind = FFTW_REDFT10;

    fftw_plan plan = fftw_plan_guru64_r2r(rank, dims, 0, …

做完FFT和IFFT后,我只能听到耳机里的噪音......这是代码:

        double* spectrum = new double[n];

        fftw_plan plan;

        plan = fftw_plan_r2r_1d(n, data, spectrum, FFTW_REDFT10, FFTW_ESTIMATE);

        fftw_execute(plan);
        fftw_destroy_plan(plan);

        plan = fftw_plan_r2r_1d(n, spectrum, data, FFTW_REDFT01, FFTW_ESTIMATE);
        fftw_execute(plan);
        fftw_destroy_plan(plan);

也许我选择了错误的FFT类型？
PS数据是初始信号

UPDATE

好的,现在代码是

        fftw_complex* spectrum = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);

        fftw_plan plan;

        plan = fftw_plan_dft_r2c_1d(n, data, spectrum, FFTW_ESTIMATE);

        fftw_execute(plan);
        fftw_destroy_plan(plan);

        plan = fftw_plan_dft_c2r_1d(n, spectrum, data, FFTW_ESTIMATE);
        fftw_execute(plan);
        fftw_destroy_plan(plan);

问题仍然存在,我的数据阵列已损坏.

更新#2

所以,问题在于我的变换大小和规范化.如果我使用真实到实际的FFTW_REDFT10和FFTW_REDFT01转换我需要使用的变换尺寸？2*N？或者是其他东西？然后我需要通过将每个元素除以2*n来规范化输出信号？
感谢大家的回复.

更新#3

感谢大家再次回复.我帮你解决了这个问题.这是工作代码:

        // FFT  
        fftw_complex* spectrum  = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * n);

        fftw_plan plan;

        plan = fftw_plan_dft_r2c_1d(n, data, spectrum, FFTW_ESTIMATE);

        fftw_execute(plan);
        fftw_destroy_plan(plan); …

我正在尝试计算FFT然后IFFT只是试试我是否可以得到相同的信号,但我不确定如何实现它.我就是这样做的FFT:

    plan = fftw_plan_r2r_1d(blockSize, datas, out, FFTW_R2HC, FFTW_ESTIMATE);
    fftw_execute(plan);

我正在尝试使用http://www.fftw.org/中的库来FFT图像,以便我可以在频域中进行卷积.但我无法弄清楚如何让它发挥作用.为了理解如何做到这一点,我试图将FFT图像转发为像素颜色数组,然后反向FFT,得到相同的像素颜色数组.这是我做的:

fftw_plan planR, planG, planB;
fftw_complex *inR, *inG, *inB, *outR, *outG, *outB, *resultR, *resultG, *resultB;

//Allocate arrays.
inR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
inB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);

outR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
outB = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);

resultR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * width * width);
resultG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * …

我正在尝试编写FORTRAN代码来评估f(r)=exp(-(r^2))使用FFTW3库的高斯函数的快速傅里叶变换.众所周知,高斯函数的傅里叶变换是另一种高斯函数.

我考虑在球面坐标中评估高斯函数的傅里叶变换积分.

因此,可以将得到的积分简化为积分[r*exp(-(r^2))*sin(kr)]dr.

我编写了以下FORTRAN代码来评估离散SINE变换DST,它是使用PURELY实数输入数组的离散傅立叶变换DFT.DST通过执行C_FFTW_RODFT00现有的FFTW3,考虑到在位置空间中的离散值是R = I*增量(i = 1,2,...,1024),以及用于DST输入数组是函数r*exp(-(r^2))NOT高斯.积分中的正弦函数是由积分[r*exp(-(r^2))*sin(kr)]dr超过SPHERICAL坐标产生的,并且它不是exp(ik.r)一般采用解析傅里叶变换时出现的虚部.

然而,结果不是动量空间中的高斯函数.

Module FFTW3
 use, intrinsic :: iso_c_binding
include 'fftw3.f03'
end module  

program sine_FFT_transform
use FFTW3
implicit none
integer, parameter :: dp=selected_real_kind(8)

real(kind=dp), parameter :: pi=acos(-1.0_dp)
integer, parameter :: n=1024 
real(kind=dp) :: delta, k
real(kind=dp) :: numerical_F_transform
integer :: i
type(C_PTR) ::  my_plan
real(C_DOUBLE), dimension(1024) :: y
real(C_DOUBLE), dimension(1024) :: yy, yk
integer(C_FFTW_R2R_KIND) :: C_FFTW_RODFT00

my_plan= fftw_plan_r2r_1d(1024,y,yy,FFTW_FORWARD, …

我有一个程序，我试图在使用 fftw 函数的 Visual Studio Code 调试器中运行该程序。它用命令编译

g++ dimer.cpp -std=c++11 -lfftw3 

在我电脑的终端上，没有抱怨未定义的引用。但是，在生成 launch.json 文件后，我的程序会抱怨 fftw 库函数和-std=c++14编译器标志。

我认为，它需要的只是额外的标志-std=c++11，并-lfftw3 在Visual Studio代码调试工作。我正在使用 Microsoft 的 C/C++ 扩展和 Code Runner 扩展。

我正在尝试将代码的 Mathematica 文档转换为 C++。

以下是我从输出中得到的错误。

Executing task: /usr/bin/g++ -g /home/msammartino/Documents/twochain/dimer.cpp -o /home/msammartino/Documents/twochain/dimer <

In file included from /usr/include/armadillo:54:0,
             from /home/msammartino/Documents/twochain/dimer.cpp:6:
/usr/include/armadillo_bits/compiler_setup.hpp:530:108: note: #pragma message: NOTE: suggest to enable C++14 mode for faster code; add -std=c++14 to compiler flags
 #pragma message ("NOTE: suggest to enable C++14 mode for faster code; …

我希望调用fftw的就地实到复杂转换函数，该函数具有以下签名：

fftw_plan fftw_plan_dft_r2c_1d(
    int n,             // transform length
    double* in,        // pointer to input array
    fftw_complex* out, // pointer to output array
    unsigned flags     // flags
);

该文档说，我应该指出我希望通过传递in和out参数的别名指针来执行就地转换。

问：怎样才能in和out不违反严格别名规则的别名？

我对特定于GCC的扩展union持开放态度（即，即使标准将其声明为未定义行为，也可以使用s进行类型校正）。即使允许此扩展，联合也不能包含动态大小的数组（在此应用程序中这是必须的-我事先不知道变换长度）。有人有什么想法吗？提前致谢。

使用FFTW(http://www.fftw.org/)而不是OpenCV的cvDFT(http://goo.gl/YCHj0)时,我能指望加速吗？

我的程序的运行时间很大程度上取决于反向和正向DFT的应用,我正在考虑使用FFTW而不是OpenCV cvDFT.

IIRC FFTW执行某种"运行时编译",而cvDFT是一个简单的FFT实现,所以我想它可以加快我的处理速度.

所以在我自己尝试之前,我想在这里问一下,如果有人偶然发现了同样的问题并解决了它.

提前致谢