FFTW使用小结(转)

最新推荐文章于 2024-10-03 17:21:09 发布
最新推荐文章于 2024-10-03 17:21:09 发布
阅读量2.6k 收藏 5
点赞数
文章标签： 算法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/Nyiragongo/article/details/129400008
版权

简介

FFTW—Fastest Fourier Transform in the West,是由 MIT 的 Matteo Frigo 博士和 Steven G. Johnson 博士开发的一个完全免费的软件包。FFTW 最初的 release 版本于 1997 年发布,最新的 release 版本 fftw-3.3.4。git路径:
https://github.com/FFTW/fftw3.git
它是一个 C 语言开发的库,支持任意大小的、任意维数的数据的离散傅里叶变换(DFT),并且还支持离散余弦变换(DCT)、离散正弦变换(DST)和离散哈特莱变换(DHT)

数据类型

FFTW 有三个版本的数据类型:double、float 和 long double，使用方法如下:

1.链接对应的库(比如 libfftw3-3、libfftw3f-3、或 ibfftw3l-3)；
2.包含同样的头文件 fftw3.h；将所有以小写"fftw_“开头的名字替换为"fftwf_”(float 版本)或"fftwl_"(long double 版本)。比如将 fftw_complex 替换为 fftwf_complex,将 fftw_execute替换为 fftwf_execute 等。
3.所有以大写"FFTW_"开头的名字不变
4.将函数参数中的 double 替换为 float 或 long double
5.最后,虽然 long double 是 C99 的标准,但你的编译器可能根本不支持该类型,或它并不能提供比 double 更高的精度。
6.fftw_malloc 考虑了数据对齐,以便使用 SIMD 指令加速,所以最好不要用 C 函数malloc 替代,而且不要将 fftw_malloc、fftw_free 和 malloc、free、 delete 等混用。
尽量使用 fftw_malloc 分配空间,而不是使用的静态数组,因为静态数组是在栈上分配的,而栈空间较小;还因为这种方式没有考虑数据对齐,不便应用SIMD 指令。
编译构造

1.默认FFTW编译生成double类型，加入参数“–enable-single”或“–enable-float”编译单精度(float)，加入参数“–enable-long-double”支持长双进度类型
2.另外,在ARM平台上可增加“–enable-neon”,使能ARM NEON流媒体加速核心
多线程

FFTW 可以多线程执行,但是多线程存在线程同步问题,这可能会降低性能。所以除非问题规模非常大,一般并不能从多线程中获益。

plan复用

用同一个 fftw_plan 执行多个数据的变换
同一个 fftw_plan 通过对输入数据赋不同值来实现不同的变换,实际上还可以利用同一个 fftw_plan 实现对不同输入输出数据的变换,也就是说可以有多个
输入输出数据数组,各自进行变换,互不影响。当然这要满足一定的条件:

输入/输出数据大小相等。
变换类型、是否原位运算不变。
对 split 数组(指实虚部分开),实部和虚部的分割方式与方案创建时相同。
数组的对齐方式相同。如果都是用 fftw_malloc 分配的则此项条件满足,除非使用FFTW_UNALIGNED 标志
总结：是用现有PLAN必须满足所有参数条件和申请PLAN时一致。

一维数据

目前使用最广泛，项目也只使用1D，并且输入数据为实数，因此变化有两种方法
1.将输入部分扩展成虚数，即所有虚数为0，然后使用一维复数变化接口

fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags);
参数说明:
n – 复数数据点数
in/out – 输入数据和输出数据，可以相同(原位相同)
sign – FFTW_FORWARD(-1)正变换 FFTW_BACKWORD(+1) 逆变换
flags – 主要有两个参数 FFTW_MEASURE/FFTW_ESTIMATE
FFTW_MEASURE 先进行预处理，对大数据(数据长度不变)连续FFT变化效果明显
FFTW_ESTIMATE 直接构造一个次最优的方案，非连续 实时选择这种方案

2.使用1D实数DFT变化，由于实数DFT变化具有Hmitian对称性，所以只需计算n/2 + 1个输出数据即可

    需要注意：
    a.如果采用原位运算则out空间必须>= 2*（n/2 + 1)
    b.1D实数DFT逆变换任何情况都会破坏IN数据
    c.当n为偶数，由Nyquust采样定理，第n/2个变化结果也是实数，所以可以把第0个和n/2个数据合并(n/2实部做0虚部)，但实际上FFTW没有这样实现，因为在多维数据时不兼容

几种方案测试

=============

使用C2C MEASURE

FFT-DFT :309us
FFT-DFT :345us
FFT-DFT :320us
FFT-DFT :349us
FFT-DFT :319us
FFT-DFT :320us
FFT-DFT :321us
FFT-DFT :324us
FFT-DFT :320us
FFT-DFT :319us
FFT-DFT :757919us
幅值:9.465116    相位:-169.040115

real    0m2.064s
user    0m2.044s
sys    0m0.020s

使用R2C

===============================
FFT-DFT :208us
FFT-DFT :195us
FFT-DFT :198us
FFT-DFT :214us
FFT-DFT :212us
FFT-DFT :212us
FFT-DFT :221us
FFT-DFT :215us
FFT-DFT :214us
FFT-DFT :216us
FFT-DFT :564592us
幅值:9.465116 相位:-169.040115

real 0m0.584s
user 0m0.552s
sys 0m0.032

使用C2C 次最优方案

===============================

复制代码
FFT-DFT :213us
FFT-DFT :204us
FFT-DFT :195us
FFT-DFT :195us
FFT-DFT :195us
FFT-DFT :196us
FFT-DFT :197us
FFT-DFT :196us
FFT-DFT :198us
FFT-DFT :197us
FFT-DFT :554788us
幅值:9.465116 相位:-169.040115

real 0m0.571s
user 0m0.544s
sys 0m0.028s
复制代码
采样点数为24000，进行1维复数/实数变化，参数选择有FFTW_MEASURE/FFTW_ESTIMATE

1.经测试，发现使用MEASURE方案时，在第一构造plan时花费约1.5秒，而一次DFT变化在300us-,构造的时间能做普通(FFTW_ESTIMATE)5000次左右，故在一般场合FFTW_ESTIMATE完全满足需求。

2.实数FFTW_ESTIMATE,时间上并没有优于复数DFT变化,只是内存开辟只有复数DFT的一半。实数DFT在构造plan时间上比复数慢，但变换比复数快，所以整体持平

3.最终选择复数和实数 ESTIMATE方案

FFTW_WISDOM

关于WISDOM,翻译自FFTW官网
FFTW -wisdom是一种实用工具来生成wisdom文件，其中包含有关如何以最佳方式计算（傅立叶）变换不同大小的保存的信息。
wisdom 的大体思路就是把生成好的策略相关的配置信息存储在磁盘里,然后在下次重新运行程序的时候,把策略相关的配置信息重新载入到内存中,这样在重新生成 plan 的时候就可以节约大量的时间。

FFTW 提供了多种方式来生成wisdom文件,使用时不必关心其中具体格式是怎么样，可以导入/出到文件/字符串等

322 FFTW_EXTERN int X(export_wisdom_to_filename)(const char *filename);
323 FFTW_EXTERN void X(export_wisdom_to_file)(FILE *output_file);
324 FFTW_EXTERN char *X(export_wisdom_to_string)(void);
325 FFTW_EXTERN void X(export_wisdom)(X(write_char_func) write_char,
326 void *data);
327 FFTW_EXTERN int X(import_system_wisdom)(void);
328 FFTW_EXTERN int X(import_wisdom_from_filename)(const char *filename);
329 FFTW_EXTERN int X(import_wisdom_from_file)(FILE *input_file);
330 FFTW_EXTERN int X(import_wisdom_from_string)(const char *input_string);
331 FFTW_EXTERN int X(import_wisdom)(X(read_char_func) read_char, void *data);

wisdom 存储起来的不是 plan 本身,而是和 plan 相关的配置信息,例如内存、寄存器等。故在导入wisdom后还是需要执行plan初始化。
example:

fftw_export_wisdom_to_filename(“wisdom.conf”);
p_str = fftw_export_wisdom_to_string();

fftw_export_wisdom_from_filename(“wisdom.conf”);
fftw_export_wisdom_from_string(p_str);

原文章链接：
https://www.cnblogs.com/zzzsj/p/14977446.html

Nyiragongo
关注
基于FFTW的快速傅里叶变换（FFT）和逆傅里叶变换（IFFT）
WunApplescript的博客
09-16 926
FFTW（Fastest Fourier Transform in the West）是一个高效的开源库，提供了各种实现FFT和IFFT的函数。通过使用FFTW库，我们可以方便地实现快速傅里叶变换和逆傅里叶变换，并在各种数字信号处理应用中应用它们。本文将介绍如何使用FFTW库实现FFT和IFFT，并给出相应的源代码和描述。需要注意的是，上述示例中的FFT和IFFT计算都采用了一维的方式。如果需要进行二维或更高维度的FFT和IFFT计算，可以使用对应的FFTW函数进行修改。，用于指定FFT的参数和操作类型。
傅里叶变换库函数FFTW
02-17
Here are FFTW lib and dll files compiled by GCC4.5 and VC2010 in Win7. But I suggest you compiling it by youself in you computer, or there will be some unpredictable errors!
FFTW介绍
qq_36314864的博客
06-09 2802
FFTW由麻省理工学院计算机科学实验室超级计算技术组开发的一套离散傅里叶变换的计算库，开源，搞笑和标准。Intel的数学库和Scililb都使用FFTW做FFT计算。1）它可计算一维或多维、实和复数据以及任意规模的2）FFTW输入数据长度任意3）FFTW支持任意多维数据4）FFTW支持SSE、SSE2 Altivec和MIPS指令集5）FFTW还包含对共享和分布式存储系统的并行变换方案：选择最佳的变换算法获得变换结果。
FFTW: 实数FFT的输出项(仅填充n/2+1项)，复数IFFT的幅度(需手工归一1/N)
最新发布
weixin_44544930的博客
10-03 543
【代码】FFTW: 实数FFT的输出项(仅填充n/2+1项)，复数IFFT的幅度(需手工归一/N)
FFTW 3 项目使用教程
gitblog_00794的博客
08-09 873
FFTW 3 项目使用教程 项目地址:https://gitcode.com/gh_mirrors/ff/fftw3 1. 项目的目录结构及介绍 FFTW 3 项目的目录结构如下： fftw3/ ├── api/ ├── doc/ ├── dft/ ├── kernel/ ├── mpi/ ├── reodft/ ├── simd-support/ ├── threads/ ├── tools/ ...
信号处理算法（4）：全球最快的傅里叶变换算法（FFTW）
雨寒sgg的博客
05-23 1万+
   本文大部分内容转载自博客【congwulong】 (https://blog.csdn.net/congwulong/article/details/7576012)   FFTW（Fastest Fourier Transform in the West）是世界上最快的FFT， 实测计算长度为10000的double数组， 单次运行时间在2ms左右。为了详细了解FFTW以及...
fftw的使用
pzs0221的专栏
05-27 4102
Qt中用fftw对音频数据变换
FFTW——一个用纯c语言写的高效FFT算法库
热门推荐
kissgoodbye2012的博客
08-01 2万+
目录 一、FFTW简介（官网翻译） 二、软件特点（官网翻译） 三、资源信息 四、VS2015下的FFTW库使用 （首先是x86版本） （然后是x64版本） 五、Qt版本下的移植 官网网址：http://www.fftw.org/ 这个已经能正常使用，不会继续研究了，这个FFTW只给预编译版本，没法学习源代码。且商用还需要花钱购买许可。 如果需要我已经移植好的源码，见如下链接：......
FFTW的使用方法 - 实数输入的数字滤波示例
weixin_38713073的博客
11-05 6818
主要记录FFTW的安装过程，FFTW的简单示例代码。并实现一个输入和输出为一维实数的数字滤波器
fftw
02-23
在提供的文件中，“fftw3.h”是FFTW库的头文件，包含了库的函数声明和数据结构定义，开发者可以通过包含这个头文件来使用FFTW的功能。"fftw3 manual.pdf"是FFTW的用户手册，详细介绍了如何编译、链接和使用FFTW库，...
FFTW介绍及FFTW库的使用
04-06
**复数类型**：FFTW使用包含两个元素的`double`型数据来表示一个复数，其中第一个元素表示实部，第二个元素表示虚部。 ```c typedef double fftw_complex[2]; ``` **精度**：FFTW默认使用双精度浮点进行运算，但用户...
Windows下FFTW库的使用
04-06
### Windows下FFTW库的使用 #### 一、FFTW简介 FFTW是Fastest Fourier Transform in the West（西方最快的傅立叶变换）的缩写，它是由麻省理工学院计算机科学实验室超级计算技术组开发的一套离散傅立叶变换...
FFTW使用手册分享.pdf
12-18
FFTW使用手册分享.pdf
FFTW.rar_FFTW
09-19
FFTW，全称为“Fastest Fourier Transform in the West”（西方最快的傅里叶变换），是一个开源的、高性能的C语言库，专门...通过熟练掌握FFTW的使用，开发者可以提升计算密集型应用的性能，尤其是在处理大规模数据时。
fftw_malloc, fftw_free 函数分析
tiankongtiankong01的博客
07-19 2922
代码分析以fftw2.15为例，原代码在fftw/malloc.c中 //空间分配函数，主要是对malloc进行了封装 void *fftw_malloc(size_t n) { char *p; int i; fftw_malloc_total += n; if (fftw_malloc_total > fftw_malloc_max) ...
【交叉编译】编译生成 x86、arm 环境下的FFTW库
challenglistic的博客
04-25 3746
【FFTW库】编译生成 x86、arm 环境下的FFTW库
【QT】采用FFTW计算FFT的频谱图和MATLAB计算不一致的情况
日常科研中的杂七杂八
12-05 1319
1、QT输出FFT频谱图2、MATLAB输出FFT频谱图。
快速傅里叶变换-FFTW库的使用-参考和翻译官方文档
KPer_Yang的博客
09-13 6748
快速傅里叶变换-FFTW库的使用-参考和翻译官方文档
使用 fopen_s读取一个wav文件并对其求WVD和STFTS并输出WVD和stft的结不使用库函数
05-20
读取 WAV 文件并解析其头部信息的代码示例如下： ```c #include <stdio.h> #include <stdlib.h> #include <stdint.h> typedef struct { char ChunkID[4]; uint32_t ChunkSize; char Format[4]; } RIFF_Chunk; typedef struct { char Subchunk1ID[4]; uint32_t Subchunk1Size; uint16_t AudioFormat; uint16_t NumChannels; uint32_t SampleRate; uint32_t ByteRate; uint16_t BlockAlign; uint16_t BitsPerSample; } FMT_Subchunk; typedef struct { char Subchunk2ID[4]; uint32_t Subchunk2Size; } DATA_Subchunk; typedef struct { RIFF_Chunk riff; FMT_Subchunk fmt; DATA_Subchunk data; } WAV_Header; WAV_Header read_wav_header(FILE* fp) { WAV_Header header; fread(&header.riff, sizeof(RIFF_Chunk), 1, fp); fread(&header.fmt, sizeof(FMT_Subchunk), 1, fp); fread(&header.data, sizeof(DATA_Subchunk), 1, fp); return header; } int main() { FILE* fp; fopen_s(&fp, "test.wav", "rb"); if (fp == NULL) { printf("Failed to open file.\n"); return -1; } WAV_Header header = read_wav_header(fp); printf("Sample rate: %u\n", header.fmt.SampleRate); printf("Bits per sample: %u\n", header.fmt.BitsPerSample); uint32_t num_samples = header.data.Subchunk2Size / (header.fmt.BitsPerSample / 8); printf("Number of samples: %u\n", num_samples); fclose(fp); return 0; } ``` 接下来，我们可以使用 FFT 算法来计算 STFT，用 Wigner-Ville 分布来计算 WVD。FFT 的代码可以使用现成的库函数（比如 FFTW），但是 WVD 的计算需要自己编写代码。以下是一个简单的 WVD 计算程序，供参考： ```c #include <stdio.h> #include <stdlib.h> #include <stdint.h> #include <math.h> #define PI 3.14159265358979323846 typedef struct { char ChunkID[4]; uint32_t ChunkSize; char Format[4]; } RIFF_Chunk; typedef struct { char Subchunk1ID[4]; uint32_t Subchunk1Size; uint16_t AudioFormat; uint16_t NumChannels; uint32_t SampleRate; uint32_t ByteRate; uint16_t BlockAlign; uint16_t BitsPerSample; } FMT_Subchunk; typedef struct { char Subchunk2ID[4]; uint32_t Subchunk2Size; } DATA_Subchunk; typedef struct { RIFF_Chunk riff; FMT_Subchunk fmt; DATA_Subchunk data; } WAV_Header; WAV_Header read_wav_header(FILE* fp) { WAV_Header header; fread(&header.riff, sizeof(RIFF_Chunk), 1, fp); fread(&header.fmt, sizeof(FMT_Subchunk), 1, fp); fread(&header.data, sizeof(DATA_Subchunk), 1, fp); return header; } double* read_wav_data(FILE* fp, uint32_t num_samples) { double* data = (double*)malloc(num_samples * sizeof(double)); if (data == NULL) { return NULL; } WAV_Header header = read_wav_header(fp); fseek(fp, sizeof(header), SEEK_SET); uint16_t bytes_per_sample = header.fmt.BitsPerSample / 8; double max_amplitude = pow(2, header.fmt.BitsPerSample - 1); for (uint32_t i = 0; i < num_samples; i++) { uint8_t buffer[4]; fread(buffer, bytes_per_sample, 1, fp); int16_t sample = 0; for (uint16_t j = 0; j < bytes_per_sample; j++) { sample |= buffer[j] << (j * 8); } data[i] = sample / max_amplitude; } return data; } double* wvd(double* data, uint32_t num_samples, uint32_t window_size, uint32_t hop_size) { double* result = (double*)malloc(num_samples * window_size * sizeof(double)); if (result == NULL) { return NULL; } double* window = (double*)malloc(window_size * sizeof(double)); if (window == NULL) { free(result); return NULL; } for (uint32_t i = 0; i < window_size; i++) { window[i] = sin(PI * (i - (window_size - 1) / 2.0) / window_size); } for (uint32_t t = 0; t < num_samples - window_size; t += hop_size) { for (uint32_t i = 0; i < window_size; i++) { for (uint32_t j = 0; j < window_size; j++) { uint32_t index = t * window_size + i * window_size + j; double real = 0.0; double imag = 0.0; for (uint32_t k = 0; k < window_size; k++) { double x = data[t + k] * window[k]; double y = data[t + i + k] * window[k]; double phase = 2 * PI * j * k / window_size; double cos_phase = cos(phase); double sin_phase = sin(phase); real += x * y * cos_phase; imag += x * y * sin_phase; } result[index] = real * real + imag * imag; } } } free(window); return result; } int main() { FILE* fp; fopen_s(&fp, "test.wav", "rb"); if (fp == NULL) { printf("Failed to open file.\n"); return -1; } WAV_Header header = read_wav_header(fp); uint32_t num_samples = header.data.Subchunk2Size / (header.fmt.BitsPerSample / 8); double* data = read_wav_data(fp, num_samples); fclose(fp); uint32_t window_size = 256; uint32_t hop_size = 128; double* WVD = wvd(data, num_samples, window_size, hop_size); free(data); free(WVD); return 0; } ``` 其中，`wvd` 函数用于计算 Wigner-Ville 分布，采用了固定长度的窗口和固定步长的重叠方式。对于每个时间点，我们计算窗口内所有时域点之间的 WVD 值，得到一个大小为 `window_size * window_size` 的矩阵，并将该矩阵展平成一个大小为 `window_size * window_size` 的向量。最终，我们得到了一个大小为 `num_samples * window_size * window_size` 的 WVD 矩阵。注意，这个程序只计算了幅度平方，如果需要计算相位信息，需要将 `real` 和 `imag` 记录下来。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Nyiragongo

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值