东方耀AI技术分享

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
热搜: 活动 交友 discuz
查看: 1224|回复: 1

[Python] 离散傅里叶变换的简单Python与C++实现

[复制链接]

1365

主题

1856

帖子

1万

积分

管理员

Rank: 10Rank: 10Rank: 10

积分
14418
QQ
发表于 2021-9-8 13:47:34 | 显示全部楼层 |阅读模式



离散傅里叶变换的简单Python与C++实现


  1. # -*- coding: utf-8 -*-
  2. __author__ = u'东方耀 微信:dfy_88888'
  3. __date__ = '2021/9/1 上午9:57'
  4. __product__ = 'PyCharm'
  5. __filename__ = '21_信号处理demo'

  6. from scipy.fftpack import fft, ifft, fftshift
  7. from scipy.signal import findfreqs, find_peaks
  8. import numpy as np
  9. import scipy
  10. from scipy.signal import get_window
  11. import matplotlib.pyplot as plt
  12. import matplotlib.ticker as mticker
  13. from matplotlib.font_manager import FontProperties

  14. font_fname = "/usr/share/fonts/wps-office/simfang.ttf"
  15. font = FontProperties(fname=font_fname)
  16. # plt.rcParams['font.sans-serif']=['simfang']#设置作图中文显示

  17. # 离散傅里叶变换的简单Python与C++实现

  18. a = np.array([1, 2, 4, 7])   # a是时域的  数字信号 采样后没有时间的概念了
  19. print("直接调包的结果:", fft(a), ifft(fft(a)))

  20. a_fft = np.zeros_like(a, dtype=np.complex)
  21. # k是频域的index   n是时域的索引index
  22. for k in range(len(a_fft)):
  23.     # a_fft[k]
  24.     for n in range(len(a)):
  25.         a_fft[k] += a[n] * np.exp(-1j * 2 * np.pi * k * n / len(a))

  26. print("a做fft的结果:", a_fft)

  27. a_fft_ifft = np.zeros_like(a_fft, dtype=np.complex)
  28. for n in range(len(a_fft_ifft)):
  29.     # a_fft_ifft[n]
  30.     for k in range(len(a_fft)):
  31.         a_fft_ifft[n] += a_fft[k] * np.exp(1j * 2 * np.pi * k * n / len(a_fft))
  32.     a_fft_ifft[n] /= len(a_fft)
  33. print("a做fft再ifft的结果:", a_fft_ifft)

  34. # 算法复杂度为 O(n**2)
复制代码


  1. #include <fftw3.h>
  2. #include <complex>
  3. #include <iostream>
  4. #include <cassert>

  5. using namespace std;

  6. typedef std::complex<float> gr_complex;   //这个float精度与fftwf一致 fftwf_complex


  7. void print_complex_array(gr_complex *p, int num){
  8.     for(int i =0; i< num;i++){
  9.         cout << "value of complex: real=" << (*(p + i)).real() << ",imag=" << (*(p + i)).imag() << endl;
  10.     }
  11.    

  12. }

  13. //FFTW编程案例:fft与ifft
  14. // 注意:target_link_libraries(main -lfftw3f)
  15. //fftw_plan一次创建可以重复使用
  16. //傅里叶变换后需要回收内存 fftw_free() 不能用 free() delete
  17. //复数的多维离散傅里叶变换  二维DFT变换  三维DFT变换

  18. int main(){
  19.     auto N = 4;
  20.     gr_complex* a;   // 数组指针
  21.     // 这个malloc还得跟数据类型对应 float
  22.     a = (gr_complex*)fftwf_malloc(sizeof(gr_complex) * N);
  23.     // 输入的数据  赋值   a = np.array([1, 2, 4, 7])
  24.     a[0] = gr_complex(1,0);
  25.     a[1] = gr_complex(2,0);
  26.     a[2] = gr_complex(4,0);
  27.     a[3] = gr_complex(7,0);
  28.     //print_complex_array(a, N);

  29.     // 变量名与Python一致
  30.     gr_complex* a_fft;
  31.     gr_complex* a_fft_ifft;
  32.     a_fft = (gr_complex*)fftwf_malloc(sizeof(gr_complex) * N);
  33.     a_fft_ifft = (gr_complex*)fftwf_malloc(sizeof(gr_complex) * N);



  34.     fftwf_plan d_fft_plan, d_ifft_plan;
  35.     cout << "C++实现离散傅里叶变换的例子[都是float类型的]!" << endl;

  36.     d_fft_plan = fftwf_plan_dft_1d(N,
  37.                                        reinterpret_cast<fftwf_complex*>(a),
  38.                                        reinterpret_cast<fftwf_complex*>(a_fft),
  39.                                        FFTW_FORWARD,
  40.                                        FFTW_ESTIMATE);
  41.     d_ifft_plan = fftwf_plan_dft_1d(N,
  42.                                     reinterpret_cast<fftwf_complex*>(a_fft),
  43.                                     reinterpret_cast<fftwf_complex*>(a_fft_ifft),
  44.                                     FFTW_BACKWARD,    // 逆变换
  45.                                     FFTW_ESTIMATE);

  46.     // 真正开始执行
  47.     fftwf_execute(d_fft_plan);
  48.     //fftw_execute(p);
  49.     cout << "fftwf变换的结果:" << endl;
  50.     print_complex_array(a_fft, N);

  51.     fftwf_execute(d_ifft_plan);
  52.     cout << "fftwf逆变换的结果:" << endl;
  53.     cout << "需注意FFTW的逆变换没有除以N,即数据正变换再反变换后是原始数据的N倍。" << endl;
  54.     print_complex_array(a_fft_ifft, N);

  55.     //收尾
  56.     fftwf_destroy_plan(d_fft_plan);
  57.     fftwf_destroy_plan(d_ifft_plan);

  58.     assert(a);
  59.     fftwf_free(a);
  60.     assert(a_fft);
  61.     fftwf_free(a_fft);
  62.     assert(a_fft_ifft);
  63.     fftwf_free(a_fft_ifft);

  64.     //assert(a == nullptr);

  65.     //先用fftwf_malloc分配输入输出内存,然后输入数据赋值,
  66.     //然后创建变换方案(fftwf_plan),然后执行变换(fftwf_execute),最后释放资源
  67.    
  68.     return 0;
  69. }
复制代码







离散傅里叶变换的公式.png
让天下人人学会人工智能!人工智能的前景一片大好!
回复

使用道具 举报

0

主题

98

帖子

200

积分

中级会员

Rank: 3Rank: 3

积分
200
发表于 2021-11-23 19:29:29 | 显示全部楼层
让天下人人学会人工智能!人工智能的前景一片大好!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|Archiver|手机版|小黑屋|人工智能工程师的摇篮 ( 湘ICP备2020019608号-1 )

GMT+8, 2024-3-29 14:04 , Processed in 0.189533 second(s), 22 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表