使用FFTW进行傅里叶变换 离散信号的一维DFT及C++编程实例

FFTW: 离散信号的一维DFT

本文介绍了如何使用FFTW信号处理库来实现离散信号的频域变换。

-不考虑相位问题。

物理含义:

 信号频率,F

 采样频率, Fs

 采样频率必须是信号频率的2倍及以上,才能保证采到的信号没有失真。

 

采样获取到数字信号后,就可以对其做DFT变换了。N个采样点,经过DFT之后,可以得到N个点的DFT结果,这N个点是以复数形式存储的。为了有利于蝶形变换运算,通常N取2的整数次方。

 

每一个点就对应着一个频率点。这个点的模值,就是该频率值下的振幅特性。具体跟原始信号的幅度有什么关系呢?

假设原始信号的峰值为A,那么DFT的结果的每个点(除了第一个点直流分量之外)的模值就是A的N/2倍。而第一个点就是直流分量,它的模值就是直流分量的N倍。

 

DFT后的N点,具有以下几个物理含义:

1' 第一个点表示0HZ,第N+1个点表示采样频率Fs(N+1个点不存在),从第1个点到N个点,这中间被N-1个点平均分成N等份,每个点的频率依次增加。例如某点n所表示的频率为:Fn=(n-1)*Fs/N

2'DFT分辨率是:Fs/N,列如,Fs=50,采样1秒钟,进行FFT,那么FFT所识别的频率是1HZ,同样是Fs=50,采样两秒钟后进行FFT,那么FFT的分辨率就是0.5HZ. 因此如果要提高频率分辨力,则必须增加采样点数,也即采样时间。频率分辨率和采样时间是倒数关系。

3' 由于采样频率是数字信号频率的两倍及以上,所以我们只需要前N/2个结果即可,从FFT的特性上来看,FFT后N个点是对称的,所以也只需要查看前N/2个结果.

https://blog.csdn.net/ssw09141324/article/details/78910939

 

FFTW信号处理库特点:

1 DFT变换后 ,幅度值为复数的平方再开方, 同时每个幅度值要除以 N/2;

2 DFT逆变换后,每个数据要除以 N;

更具体可以看参考的博客,这里就不重复了。

参考:

https://blog.csdn.net/ssw09141324/article/details/78910939

https://blog.csdn.net/grafx/article/details/38750107

 

C++代码:

    int N = 8;
    double* in = NULL;
    // 如果要使用float版本,需先引用float版本的lib库,然后在fftw后面加上f后缀即可.
    fftw_complex* out = NULL;// fftwf_complex --> 即为float版本
    fftw_plan p;
    in = (double*)fftw_malloc(sizeof(double) * N);
    out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * N);

    double F_A = 1.0 / N;//信号A频率
    double F_B = 1.0 / N * 2;//信号B频率
    double Fs = 2;//hz 采样频率;
    double timeVal = 0;

    // 输入纯实数
    printf("初始信号-时域:");
    for (int i = 0; i < N; i++)
    {
        timeVal = (1.0 / Fs) * i;
        in[i] = sin(2 * PI * F_A * timeVal) + sin(2 * PI * F_B * timeVal);
        printf("%.2f ", in[i]);
    }
    printf("\n\n");

    // 傅里叶变换
    p = FFTW3_H::fftw_plan_dft_r2c_1d(N, in, out, FFTW_ESTIMATE);
    fftw_execute(p);
    
    double dx3 = Fs / N;//频率轴最小刻度单位
    // 输出幅度谱
    printf("频域(频率 - 幅度):");
    for (int i = 0; i < N; i++)
    {
        double val = sqrt(out[i][0] * out[i][0] + out[i][1] * out[i][1]);
        printf("(%.3f - %.2f) ", dx3 * i, val / (N / 2));
    }
    printf("\n\n");

    for (size_t i = 0; i < N; i++)
    {
        in[i] = 0;
    }

    //逆变换
    p = FFTW3_H::fftw_plan_dft_c2r_1d(N, out, in, FFTW_ESTIMATE);
    fftw_execute(p);

    //信号还原
    printf("逆变换 时域:");
    for (int i = 0; i < N; i++)
    {
        printf("%.2f ", in[i] / N);
    }
    printf("\n");

    // 释放资源
    fftw_destroy_plan(p);
    fftw_free(in);
    fftw_free(out);

 

发布了7 篇原创文章 · 获赞 4 · 访问量 3321
展开阅读全文

没有更多推荐了,返回首页

©️2019 CSDN 皮肤主题: 编程工作室 设计师: CSDN官方博客

分享到微信朋友圈

×

扫一扫,手机浏览