Jos*_*fer 30
我使用Sliding DFT,在每次样本到达输入缓冲区时需要进行傅里叶变换的情况下,它比FFT快许多倍.
这是基于以下事实:一旦您对最后N个样本执行了傅里叶变换,并且新样本到达,您可以"撤消"最旧样本的效果,并在一次通过中应用最新样本的效果通过傅里叶数据!这意味着对于FFT ,滑动DFT性能是O(n)与O(Log2(n)乘以n)相比.此外,对于保持性能的缓冲区大小,对2的幂没有限制.
下面的完整测试程序将滑动DFT与fftw进行比较.在我的生产代码中,我将以下代码优化为不可读性,使其速度提高了三倍.
#include <complex>
#include <iostream>
#include <time.h>
#include <math_defines.h>
#include <float.h>
#define DO_FFTW // libfftw
#define DO_SDFT
#if defined(DO_FFTW)
#pragma comment( lib, "d:\\projects\\common\\fftw\\libfftw3-3.lib" )
namespace fftw {
#include <fftw/fftw3.h>
}
fftw::fftw_plan plan_fwd;
fftw::fftw_plan plan_inv;
#endif
typedef std::complex<double> complex;
// Buffer size, make it a power of two if you want to improve fftw
const int N = 750;
// input signal
complex in[N];
// frequencies of input signal after ft
// Size increased by one because the optimized sdft code writes data to freqs[N]
complex freqs[N+1];
// output signal after inverse ft of freqs
complex out1[N];
complex out2[N];
// forward coeffs -2 PI e^iw -- normalized (divided by N)
complex coeffs[N];
// inverse coeffs 2 PI e^iw
complex icoeffs[N];
// global index for input and output signals
int idx;
// these are just there to optimize (get rid of index lookups in sdft)
complex oldest_data, newest_data;
//initilaize e-to-the-i-thetas for theta = 0..2PI in increments of 1/N
void init_coeffs()
{
for (int i = 0; i < N; ++i) {
double a = -2.0 * PI * i / double(N);
coeffs[i] = complex(cos(a)/* / N */, sin(a) /* / N */);
}
for (int i = 0; i < N; ++i) {
double a = 2.0 * PI * i / double(N);
icoeffs[i] = complex(cos(a),sin(a));
}
}
// initialize all data buffers
void init()
{
// clear data
for (int i = 0; i < N; ++i)
in[i] = 0;
// seed rand()
srand(857);
init_coeffs();
oldest_data = newest_data = 0.0;
idx = 0;
}
// simulating adding data to circular buffer
void add_data()
{
oldest_data = in[idx];
newest_data = in[idx] = complex(rand() / double(N));
}
// sliding dft
void sdft()
{
complex delta = newest_data - oldest_data;
int ci = 0;
for (int i = 0; i < N; ++i) {
freqs[i] += delta * coeffs[ci];
if ((ci += idx) >= N)
ci -= N;
}
}
// sliding inverse dft
void isdft()
{
complex delta = newest_data - oldest_data;
int ci = 0;
for (int i = 0; i < N; ++i) {
freqs[i] += delta * icoeffs[ci];
if ((ci += idx) >= N)
ci -= N;
}
}
// "textbook" slow dft, nested loops, O(N*N)
void ft()
{
for (int i = 0; i < N; ++i) {
freqs[i] = 0.0;
for (int j = 0; j < N; ++j) {
double a = -2.0 * PI * i * j / double(N);
freqs[i] += in[j] * complex(cos(a),sin(a));
}
}
}
double mag(complex& c)
{
return sqrt(c.real() * c.real() + c.imag() * c.imag());
}
void powr_spectrum(double *powr)
{
for (int i = 0; i < N/2; ++i) {
powr[i] = mag(freqs[i]);
}
}
int main(int argc, char *argv[])
{
const int NSAMPS = N*10;
clock_t start, finish;
#if defined(DO_SDFT)
// ------------------------------ SDFT ---------------------------------------------
init();
start = clock();
for (int i = 0; i < NSAMPS; ++i) {
add_data();
sdft();
// Mess about with freqs[] here
//isdft();
if (++idx == N) idx = 0; // bump global index
if ((i % 1000) == 0)
std::cerr << i << " iters..." << '\r';
}
finish = clock();
std::cout << "SDFT: " << NSAMPS / ((finish-start) / (double)CLOCKS_PER_SEC) << " fts per second." << std::endl;
double powr1[N/2];
powr_spectrum(powr1);
#endif
#if defined(DO_FFTW)
// ------------------------------ FFTW ---------------------------------------------
plan_fwd = fftw::fftw_plan_dft_1d(N, (fftw::fftw_complex *)in, (fftw::fftw_complex *)freqs, FFTW_FORWARD, FFTW_MEASURE);
plan_inv = fftw::fftw_plan_dft_1d(N, (fftw::fftw_complex *)freqs, (fftw::fftw_complex *)out2, FFTW_BACKWARD, FFTW_MEASURE);
init();
start = clock();
for (int i = 0; i < NSAMPS; ++i) {
add_data();
fftw::fftw_execute(plan_fwd);
// mess about with freqs here
//fftw::fftw_execute(plan_inv);
if (++idx == N) idx = 0; // bump global index
if ((i % 1000) == 0)
std::cerr << i << " iters..." << '\r';
}
// normalize fftw's output
for (int j = 0; j < N; ++j)
out2[j] /= N;
finish = clock();
std::cout << "FFTW: " << NSAMPS / ((finish-start) / (double)CLOCKS_PER_SEC) << " fts per second." << std::endl;
fftw::fftw_destroy_plan(plan_fwd);
fftw::fftw_destroy_plan(plan_inv);
double powr2[N/2];
powr_spectrum(powr2);
#endif
#if defined(DO_SDFT) && defined(DO_FFTW)
// ------------------------------ ---------------------------------------------
const double MAX_PERMISSIBLE_DIFF = 1e-11; // DBL_EPSILON;
double diff;
// check my ft gives same power spectrum as FFTW
for (int i = 0; i < N/2; ++i)
if ( (diff = abs(powr1[i] - powr2[i])) > MAX_PERMISSIBLE_DIFF)
printf("Values differ by more than %g at index %d. Diff = %g\n", MAX_PERMISSIBLE_DIFF, i, diff);
#endif
return 0;
}
Sri*_*ram 16
如果在一个图中需要幅度,频率和时间,则变换称为时频分解.最受欢迎的称为短时傅里叶变换.它的工作原理如下:
1.取一小部分信号(比如1秒)
2.用一个小窗口(比如说5毫秒)
对其进行窗口化3.计算窗口信号的一维傅立叶变换.
4.将窗口移动少量(2.5 ms)
5.重复上述步骤直到信号结束.
6.所有这些数据都输入一个矩阵,然后用于创建信号的3D表示,显示其沿频率,幅度和时间的分解.
窗口的长度将决定您在频域和时域中可以获得的分辨率.请点击这里关于STFT更多详细信息,搜索"ROBI Polikar"的小波变换的教程对于一个外行的介绍上面.
编辑1:
你采用了一个窗口函数(那里有无数的函数 - 这是一个列表.最直观的是一个矩形窗口,但最常用的是汉明/汉宁窗口函数.如果你有一个,你可以按照下面的步骤手中的纸笔和绘制它.
假设您获得的信号长1秒并命名x[n].窗口函数长5毫秒,并命名w[n].将窗口放在信号的开头(因此窗口的末端与信号的5ms点重合)并将其相乘x[n],w[n]如下所示:
y[n] = x[n] * w[n]- 信号的逐点乘法.
进行FFT计算y[n].
然后你将窗口移动少量(比如2.5毫秒).所以现在窗口从2.5ms延伸到7.5ms的信号x[n].重复乘法和FFT生成步骤.换句话说,您有2.5毫秒的重叠.您将看到更改窗口长度和重叠会在时间轴和频率轴上提供不同的分辨率.
完成此操作后,您需要将所有数据提供到矩阵中,然后显示它.重叠用于最小化在边界处可能出现的误差,并且还在这样的短时间帧内获得更一致的测量.
PS:如果你已经理解了信号的STFT和其他时频分解,那么你应该对第2步和第4步没有任何问题.你还没有理解上面提到的步骤让我觉得你应该重新考虑时频分解也.