QT下利用FFTW实现功率谱密度分析（周期图法）

写在前面

出于项目需求，需要在QT的msvc 2015编译环境下实现功率谱密度估计。

所以这个需求很直接的关键词就是

1. C++

2. QT

3. 傅里叶变换

4. 功率谱密度估计

但是在网上搜索了很久，没有一个非常直接的答案。所以我觉得很有必要写一篇文章来说明，在QT msvc2015 32bit环境下如何实现功率谱密度估计。

本文并没有经过严格编辑，但代码是已经可以跑通的，如果我的说明描述存在问题，请见谅。

环境 or 需求

Win7 64bit

QT msvc2015 32bit

FFTW3 ( fftw-3.3.5-dll32 )  官网：http://www.fftw.org/

QCustomPlot 官网：https://www.qcustomplot.com

正式开始

fftw3、qcustomplot、qt的使用方法以及功率谱密度估计的原理在这里就不赘述了，我们直奔主题——如何用代码实现。

原理

对于一个信号的功率谱密度估计，我采用的是originlab当中使用的原理，如下图：

(链接：https://www.originlab.com/doc/Origin-Help/FFT1-Algorithm)

代码

#include "mainwindow.h"
#include "ui_mainwindow.h"
#include "fftw3.h"
#include <QDebug>

MainWindow::MainWindow(QWidget *parent) :
    QMainWindow(parent),
    ui(new Ui::MainWindow)
{
    ui->setupUi(this);
    int N = 10000; //采集点数目
    int Fs = 100000; //采集频率，这里是100kHz，也就是每隔10us采集一个点
    double sens = 1.90648; //这里是我项目里面的一个灵敏度因子，可以当作一个常数来看待

    QString fileName = "C:/Users/CMCS-LZC/Desktop/QT_vs_test_FFT/FFT_test/cmpms_test_D.txt"; //选择路径
    QFile file(fileName);
    QVector<double> array; //用来保存信号
    QVector<double> array_x; //画图的时候使用，array_x作为横坐标，array作为纵坐标
    if(file.open(QIODevice::ReadOnly))
    {
        QTextStream stream(&file);
        while(!file.atEnd())
        {
            double buf;
            QStringList list=stream.readAll().split("\n");
            QListIterator<QString> li(list);
            while(li.hasNext())
            {
                buf=li.next().toDouble();
                double temp = (buf/sens) * 1e-6 ;
                array.append(temp);
            }
        }
    }

    //进行傅里叶变换
    fftw_complex *in, *out;
    in = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * N);

    for (int i = 0; i < N; i++)
    {
        in[i][0] = array[i];
        in[i][1] = 0;
        array_x.append(i);
    }

    fftw_plan my_plan;
    my_plan = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(my_plan);

    QVector<double> data_power;
    QVector<double> data_power_x;
    QVector<double> data_freq;

    //功率谱密度估计
    for( int i = 82; i < 4905; i++) //这里的i的区间可以自由选择。
    {
        data_power.append((out[i][0] * out[i][0] + out[i][1] * out[i][1])/(Fs*N));
//        data_cali_psd_frequency.append((2*Pi*i)/(N*T));
        data_freq.append(i*10);
        data_power_x.append(i-82);
    }
    fftw_destroy_plan(my_plan);
    fftw_free(in);
    fftw_free(out);
    //傅里叶变换结束，功率谱密度估计计算结束

    // 求积分来计算弹簧系数
    int length_psd = data_power.length();
    double area = 0;
    for(int i = 0; i<(length_psd-1); i++)
    {
        area = area + 0.5 *( (data_freq[i+1]-data_freq[i]) * (data_power[i] + data_power[i+1]) );
    }
    qDebug()<<"area is "<<area;

    // 计算弹簧系数
    double k_b = 1.38 * 1e-23;
    double Temperature = 300;
    double springconstant = 0.764 * (k_b * Temperature ) / area ; // N/m ，锥形针0.764，矩形针0.817
    qDebug()<<"Spring Constant is "<<springconstant*1e3 <<"pN/nm";

    ui->plot_1->addGraph();
    ui->plot_1->graph(0)->setData(array_x, array);
    ui->plot_1->rescaleAxes(true);

    ui->plot_2->addGraph();
    ui->plot_2->graph(0)->setData(data_freq, data_power);
    ui->plot_2->rescaleAxes(true);

}

MainWindow::~MainWindow()
{
    delete ui;
}

运行结果

 结语

本文所做的功率谱密度估计使用的是周期图法，是实现起来较为简单的一种方法。

有任何问题，请说明来意，email（lzc2079@qq.com）给我。