傅里叶变换

正弦特征

正余弦拥有其他信号所不具备的性质：正弦曲线具有保真度（正弦曲线信号输入后，输出的仍是正弦曲线，只有幅度和相位可能发生变化，但是频率和波的形状仍是一样的）

原理

任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。

时域信号在经过傅立叶变换的分解之后，变为了多个不同正弦波信号的叠加，我们再去分析这些正弦波的频率。这样就成功将信号时域分析切换到频域分析。

对于离散信号（在连续信号上采样得到的信号）的变换只有离散傅立叶变换（DFT）才能被适用，且对于计算机来说只有离散的和有限长度的数据才能被处理。同时计算DFT的复杂度是N^2，而进行FFT（快速傅里叶变换）的运算复杂度是N*log(N)。----> 编程实现的傅里叶变换实际是FFT。

基 2 FFT算法

离散傅里叶变换（DFT）


参考链接: FFT变换-北邮张建华老师







傅里叶变换的核心：对输入信号持续进行奇偶二分

FFT总体思路

输入X(n)（时域数据）–>> 序列倒序 -->> 蝶形运算 -->> 输出结果
注：AD采集数据的点数必须为2^n；

编程实现

参考链接: 快速傅里叶变换
Complex.h文件

#ifndef COMPLEX_H
#define COMPLEX_H

#include <math.h>

//声明复数结构体
struct Complex
{
	double real;//实部
	double imag;//虚部
}

//复数相加
Complex Complex_Add(Complex a, Complex b)
{
	Complex ret;
	ret.real = a.real + b.real;
	ret.imag = a.imag + b.imag;
	return ret;
}

//复数相减
Complex Complex_sub(Complex a, Complex b)
{
	Complex ret;
	ret.real = a.real - b.real;
	ret.imag = a.imag - b.imag;
	return ret;
}

//复数相乘
Complex Complex_Mult(Complex a, Complex b)
{
	Complex ret;
	ret.real = a.real*b.real - a.imag*b.imag;
	ret.imag = a.imag*b.real + a.real*b.imag;
	return ret;
}

//复数求模
double Complex_Model(Complex a)
{
	return sqrt(pow(a.real, 2) + pow(a.imag, 2));
}

//将时域数据转换为复数形式表示
//sample:数据转换为复数形式保存
//buf:时域原始数据
//num:原始数据个数
void Complex_Set(Complex *sample, double *buf, int num)
{
	int i;
	for(i = 0;i < num;i++)
	{
		sample[i].real = buf[i];
		sample[i].imag = 0;
	}
}

#endif

FFT.cpp文件

#include <stdio.h>
#include <stdlib.h>
#include "Complex.h"

double PI = 4.0*atan(1);//定义pi，tan(pi/4)=1增加pi的精度
//举例数组
double data[] = {1,2,3,4,5,6,7,8};

void FFTButterfly(double *source, int len)
{
	int i,j,k,l,p;
	int index;//数组下标
	int gap;//间隔
	int increment;//增量
	Complex factor;//旋转因子
	Complex tmp1;//蝶形运算中暂存
	//声明时域数据复数形式存储空间
	Complex *TimeD = (Complex *)malloc(sizeof(Complex)*len);
	//声明序列倒序存储空间 + FFT运算结果存储
	Complex *TimeDInvert = (Complex *)malloc(sizeof(Complex)*len);
	
	//时域数据转为复数形式
	Complex_Set(TimeD, data, len);
	
	//求len点的FFT分级数
	int level = log2(len);
	
	//序列倒序
	for(i = 0; i < len; i++)
	{
		k = 0;
		for(j = 0; j < level; j++)
		{
			if(((i >> j) & 1) == 1)
			{
				k += 1 << (level - j - 1);
			}
		}
		TimeDInvert[k] = TimeD[i];
	}

	//蝶形运算
	for(l = 1;l <= level;l++)
	{
		int gap = int(pow(2, (l - 1)));//间隔
		int increment = int(pow(2, (level - l)));//增量
		for(i = 0;i < increment;i++)//蝶形运算的次数
		{
			for(j = 0;j < gap;j++)//每个蝶形运算的乘法次数
			{
				index = j + 2 * gap * i;//数组下标
				p = j * increment;
				factor.real = cos(2 * PI * p / double(len));
				factor.imag = -sin(2 * PI * p / double(n));
				tmp1 = Complex_Mult(factor, TimeDInvert[index + gap]);
				TimeDInvert[index + gap] = Complex_sub(TimeDInvert[index], tmp1);
				TimeDInvert[index] = Complex_Add(TimeDInvert[index], tmp1);
			}
		}
	}

	//打印结果
	for(i = 0;i < len;i++)
	{
		printf("%.4f ",Complex_Model(TimeDInvert[i]));
		printf("%.4f + ",TimeDInvert[i].real);
		printf("%.4f j\n",TimeDInvert[i].imag);
	}

	free(TimeD);
	free(TimeDInvert);
}

int main()
{
	FFTButterfly(data, sizeof(data)/sizeof(double));
	system("pause");
	return 0;
}

Matlab运算结果

窗函数

FFT变换是在完美情况下（即认为被处理的信号为周期信号），但在实际情况下AD采集的数据是离散的，对离散数据进行周期扩展时，会出现数据跳变的情况（数据跳变 <–> 相位不连续，相当于引入相位调制）生成额外的频率成分，这种现象称为频谱泄露
窗函数解释: 参考链接

窗函数C语言实现: 参考链接

使用方法

1）时间窗多项式的样点数与AD采集的数据样本数相同。
2）将时间窗的系数与存放时域数值的数组元素一一相乘。

分帧

链接: 分帧原理
链接: C代码实现