回顾
上一节我们已经讲过了窗函数,现在信号已经看起来像是一个周期函数了,通过合理选择窗函数,我们既成功的保证了数据精度,又减少了频谱泄露!
现在,信号要经过FFT,数据取半和求模啦!
开始
第一点注意
输入FFT的点数必须是2的整数幂次。
也就是像1024,2048,4096这样的点。
可是,如果我们采样的点数就不是2的整数次幂怎么办呢?末尾补零!
举个例子:
我们有一串信号有4000个点,那么就要在后面补96个零再过FFT。
下面我们进一段程序:
uint32_t WaveMeasureUtils::findNextPow2(uint32_t v)
{
v--;
v |= v >> 1;
v |= v >> 2;
v |= v >> 4;
v |= v >> 8;
v |= v >> 16;
v++;
return v;
}
使用这段找下一幂次的程序程序,输入数据点的个数,它能够帮你找到需要进行的FFT点数,我们也就能求得补零的个数了。
经过了补零,我们可以愉快的经过FFT啦!再上一段工作良好的FFT的程序,给大家做参考(使用Qt编写,需要按情况进行改写):
.h:
#ifndef FFT_H
#define FFT_H
#include <QObject>
#include <QDebug>
/*
*
*
* 计算傅里叶变换频谱
*
*
* */
#define MAX_MATRIX_SIZE 4194304 // 2048 * 2048
#define PI 3.141592653
#define MAX_VECTOR_LENGTH 10000
typedef struct Complex
{
double rl; //实部
double im; //虚部
}Complex;
class fft : public QObject
{
Q_OBJECT
public:
explicit fft(QObject *parent = nullptr);
//傅里叶转换 频域
bool fft1(QVector<Complex>inVec, const int len, QVector<Complex> &outVec);
int get_computation_layers(int num);
bool is_power_of_two(int num);
void test();
signals:
public slots:
};
#endif // FFT_H
.cpp:
#include "fft.h"
#include <QDebug>
fft::fft(QObject *parent) : QObject(parent)
{
}
bool fft::fft1(QVector<Complex> inVec, const int len, QVector<Complex> &outVec)
{
if ((len <= 0) || (inVec.isEmpty()) || ( outVec.isEmpty()))
return false;
if (!is_power_of_two(len))
return false;
Complex *pVec = new Complex[len];
Complex *Weights = new Complex[len];
Complex *X = new Complex[len];
int *pnInvBits = new int[len];
//memcpy(pVec, inVec, len*sizeof(Complex));
for(int i = 0; i < len;i++)
{
pVec[i].im = inVec.at(i).im;
pVec[i].rl = inVec.at(i).rl;
}
// 计算权重序列
double fixed_factor = (-2 * PI) / len;
for (int i = 0; i < len / 2; i++) {
double angle = i * fixed_factor;
Weights[i].rl = cos(angle);
Weights[i].im = sin(angle);
}
for (int i = len / 2; i < len; i++) {
Weights[i].rl = -(Weights[i - len / 2].rl);
Weights[i].im = -(Weights[i - len / 2].im);
}
int r = get_computation_layers(len);
// 计算倒序位码
int index = 0;
for (int i = 0; i < len; i++) {
index = 0;
for (int m = r - 1; m >= 0; m--) {
index += (1 && (i & (1 << m))) << (r - m - 1);
}
pnInvBits[i] = index;
X[i].rl = pVec[pnInvBits[i]].rl;
X[i].im = pVec[pnInvBits[i]].im;
}
// 计算快速傅里叶变换
for (int L = 1; L <= r; L++) {
int distance = 1 << (L - 1);
int W = 1 << (r - L);
int B = len >> L;
int N = len / B;
for (int b = 0; b < B; b++) {
int mid = b*N;
for (int n = 0; n < N / 2; n++) {
int index = n + mid;
int dist = index + distance;
pVec[index].rl = X[index].rl + (Weights[n*W].rl * X[dist].rl - Weights[n*W].im * X[dist].im); // Fe + W*Fo
pVec[index].im = X[index].im + Weights[n*W].im * X[dist].rl + Weights[n*W].rl * X[dist].im;
}
for (int n = N / 2; n < N; n++) {
int index = n + mid;
pVec[index].rl = X[index - distance].rl + Weights[n*W].rl * X[index].rl - Weights[n*W].im * X[index].im; // Fe - W*Fo
pVec[index].im = X[index - distance].im + Weights[n*W].im * X[index].rl + Weights[n*W].rl * X[index].im;
}
}
for(int i = 0; i< len;i++)
{
X[i].im = pVec[i].im;
X[i].rl = pVec[i].rl;
}
}
for(int i = 0; i < len;i++)
{
outVec[i].im = pVec[i].im;
outVec[i].rl = pVec[i].rl;
}
if (Weights) delete[] Weights;
if (X) delete[] X;
if (pnInvBits) delete[] pnInvBits;
if (pVec) delete[] pVec;
return true;
}
int fft::get_computation_layers(int num)
{
int nLayers = 0;
int len = num;
if (len == 2)
return 1;
while (true) {
len = len / 2;
nLayers++;
if (len == 2)
return nLayers + 1;
if (len < 1)
return -1;
}
}
bool fft::is_power_of_two(int num)
{
int temp = num;
int mod = 0;
int result = 0;
if (num < 2)
return false;
if (num == 2)
return true;
while (temp > 1)
{
result = temp / 2;
mod = temp % 2;
if (mod)
return false;
if (2 == result)
return true;
temp = result;
}
return false;
}
第二点注意
FFT的输出数据是要取半的!
过了FFT,我们得到了同样多个数的数据点。
但是!我们只能取他们的一半。根据采样定理,频率高于0.5采样率的数据点过FFT是无法解析正常的。
FFT的数据点从第0个到最后一个一一对应着0~采样率,频率分辨率为:采样率/FFT点数
举个例子:如果你的采样率是20khz,点数是4096,那么在第2047点的时候就是10khz的频率对应的点了,再往后的点数没有可用价值。
实际上,如果非要画出来,数据点是中心对称的:
第三点注意
FFT的输出是实实在在的复数
FFT输入的点是个实数(即使输入要求填写虚部,那虚部就为零)
而FFT的输出是实实在在的复数(实部虚部均有用处)
我们要求信号的各种幅值的时候,就需要对他们进行取模运算(勾股定理)。
for(auto i = 0;i<out_.size()/2;i++)
{
x1<<i;
//进行取模
y1<<sqrt(out_[i].rl * out_[i].rl+out_[i].im * out_[i].im);
}
到这里,FFT和数据取半,求模就结束了,本文图片来源:https://www.jianshu.com/p/b2f64234cc3b,在此声明。
下一篇我们聊聊乘恢复系数和除不同系数。
欢迎大家关注,点赞哦,您的支持是我加速更新和为您写出更好文章的动力!