3.BCI中的信号处理——快速傅里叶变换

WR188

已于 2024-04-03 21:11:10 修改

阅读量581

点赞数 16

分类专栏：脑机文章标签：信号处理人机交互

于 2024-04-01 15:31:22 首次发布

本文链接：https://blog.csdn.net/weixin_45041480/article/details/137234085

版权

脑机专栏收录该内容

3 篇文章 0 订阅

订阅专栏

傅里叶变换（FT）的基本思想是将一个定义在区间[-T/2,T/2]之间的时变信号 $s(t)$ 分解为不同频率 $\omega$ 的正弦和余弦波的加权和：

$s(t)=\frac{a_0}{2}+\sum_{n=1}^{\infty }a_n\cos(n\omega t)+\sum_{n=1}^{\infty }b_n\sin(n\omega t)$

在实际工程应用中，通常是在离散时间间隔上对大脑信号进行采样，因此需要对连续时间傅里叶变换进行修改，以满足离散采样信号的需求。对于输入时间序列 $S(t)$ ，仅在时间点 $t=0,\cdots,T-1$ 有采样值，离散傅里叶变换（DFT）将该序列转换成相应复数形式的傅里叶系数：

$C(n)=\frac{1}{T}\sum_{t=0}^{T-1}S(t)e^{-jn\omega t},n=1,\cdots,T-1$

DFT对于 $T$ 个数据的计算复杂度是 $T^2$ ，当 $T$ 较大时，运算会十分缓慢。本部分引入快速傅里叶变换（FFT），FFT的计算复杂度是 $T\log(T)$ 。FFT数值计算方法参考博客，本部分提供一个基于C/C++的FFT实例，头文件“FFT.h”可作为工具直接使用，针对离散时序信号 $S(t)=\sin(0.1t)$ 分析其频域特征。除此之外，fftw库提供了经过优化的FFT计算工具。

FFT.h：

//本文件的目的是实现一个快速傅里叶变换（FFT）的功能
#ifndef FFT_H
#define FFT_H
#include<iostream>
#include<cmath>
#include<complex>
#include<vector>
class FFT{
    public:
    //std::complex<double> var{a,b}是构造复数（a+bi），支持拷贝构造
    static std::complex<double> cis(double angle){//输入角度为弧度制
        std::complex<double> var{std::cos(angle),std::sin(angle)};
        return var;
    }
    static std::vector<std::complex<double>> fft(std::vector<double> input){
        double pi=3.14159265359; 
        int shape=input.size();
        if(shape<2) return {std::complex<double>(input[0],0)};
        std::vector<double> seq0,seq1;
        for(int i=0;i<shape;i++){
            if(i%2==0) seq0.push_back(input[i]);
            else seq1.push_back(input[i]);
        }
        std::vector<std::complex<double>> seq0_fft = fft(seq0);
        std::vector<std::complex<double>> seq1_fft = fft(seq1);

        std::vector<std::complex<double>> result(shape);
        for (int k =0; k<shape/2; ++k) {
            double tk=k*2*pi/shape;
            result[k] = seq0_fft[k]+seq1_fft[k]*cis(-tk);
            result[k+shape/2]=seq0_fft[k]-seq1_fft[k]*cis(-tk);
        }
        return result;
    }
};

#endif

main.cpp:

#include <iostream>
#include <vector>
#include <complex>
#include "FFT.h"
int main(){
    FFT fft_test;
    std::vector<double> input;
    for(int i=0;i<1000;i++){
        input.push_back(std::sin(i*0.1));
    }
   std::vector<std::complex<double>> output = fft_test.fft(input);
   std::cout<<"print FFT result:"<<std::endl;
   for(unsigned j=0;j<output.size();j++)    std::cout<<output[j]<<std::endl;
    return 0;
}