ifft 结果与原始信号不同
ifft results are different from original signal
FFT 工作正常,但当我想使用 IFFT 时,我总是从其结果中看到相同的图形。结果很复杂,无论原始信号如何,图形始终相同。
实部图是一个 -sin,周期 = 帧大小
虚部是同周期的-cos
哪里会出问题?
原始信号:
IFFT真实值(图片只有一半帧):
我用的算法FFT。
double** FFT(double** f, int s, bool inverse) {
if (s == 1) return f;
int sH = s / 2;
double** fOdd = new double*[sH];
double** fEven = new double*[sH];
for (int i = 0; i < sH; i++) {
int j = 2 * i;
fOdd[i] = f[j];
fEven[i] = f[j + 1];
}
double** sOdd = FFT(fOdd, sH, inverse);
double** sEven = FFT(fEven, sH, inverse);
double**spectr = new double*[s];
double arg = inverse ? DoublePI / s : -DoublePI / s;
double*oBase = new double[2]{ cos(arg),sin(arg) };
double*o = new double[2]{ 1,0 };
for (int i = 0; i < sH; i++) {
double* sO1 = Mul(o, sOdd[i]);
spectr[i] = Sum(sEven[i], sO1);
spectr[i + sH] = Dif(sEven[i], sO1);
o = Mul(o, oBase);
}
return spectr;
}
"butterfly" 部分错误地应用了系数:
for (int i = 0; i < sH; i++) {
double* sO1 = sOdd[i];
double* sE1 = Mul(o, sEven[i]);
spectr[i] = Sum(sO1, sE1);
spectr[i + sH] = Dif(sO1, sE1);
o = Mul(o, oBase);
}
旁注:
我保留了你的注释,但这让事情变得混乱:
fOdd
有索引 0, 2, 4, 6, ... 所以它应该是 fEven
fEven
有索引 1, 3, 5, 7, ... 所以它应该是 fOdd
实际上 sOdd
应该是 sLower
而 sEven
应该是 sUpper
因为它们对应于 0:s/2
和 s/2:s-1
元素光谱分别为:
sLower = FFT(fEven, sH, inverse); // fEven is 0, 2, 4, ...
sUpper = FFT(fOdd, sH, inverse); // fOdd is 1, 3, 5, ...
那么蝴蝶就变成了:
for (int i = 0; i < sH; i++) {
double* sL1 = sLower[i];
double* sU1 = Mul(o, sUpper[i]);
spectr[i] = Sum(sL1, sU1);
spectr[i + sH] = Dif(sL1, sU1);
o = Mul(o, oBase);
}
这样写比较容易与 pseudocode example on wikipedia.
进行比较
@Dai 是正确的,你会泄漏大量内存
关于内存,您可以使用 std::vector
来封装动态分配的数组,并确保在执行离开范围时释放它们。您可以使用 unique_ptr<double[]>
但性能提升在 IMO 上不值得,并且您失去了 at()
方法的安全性。
(基于@Robb 的回答)
其他一些提示:
- 避免神秘的标识符 - 程序应该是可读的,像“
f
”和“s
”这样的名称会使您的程序更难阅读和维护。
- 基于类型的匈牙利符号不受欢迎,因为现代编辑器会自动显示类型信息,因此它给标识符名称增加了不必要的复杂性。
- 索引使用
size_t
,而不是int
- STL 是你的朋友,使用它!
- 通过使用
const
预防只读数据的意外突变来预防错误。
像这样:
#include <vector>
using namespace std;
vector<double> fastFourierTransform(const vector<double> signal, const bool inverse) {
if( signal.size() < 2 ) return signal;
const size_t half = signal.size() / 2;
vector<double> lower; lower.reserve( half );
vector<double> upper; upper.reserve( half );
bool isEven = true;
for( size_t i = 0; i < signal.size(); i++ ) {
if( isEven ) lower.push_back( signal.at( i ) );
else upper.push_back( signal.at( i ) );
isEven = !isEven;
}
vector<double> lowerFft = fastFourierTransform( lower, inverse );
vector<double> upperFft = fastFourierTransform( upper, inverse );
vector<double> result;
result.reserve( signal.size() );
double arg = ( inverse ? 1 : -1 ) * ( DoublePI / signal.size() );
// Ideally these should be local `double` values passed directly into `Mul`.
unique_ptr<double[]> oBase = make_unique<double[]>( 2 );
oBase[0] = cos(arg);
oBase[1] = sin(arg);
unique_ptr<double[]> o = make_unique<double[]>( 2 );
o[0] = 0;
o[1] = 0;
for( size_t i = 0; i < half; i++ ) {
double* lower1 = lower.at( i );
double* upper1 = Mul( o, upper.at( i ) );
result.at( i ) = Sum( lower1, upper1 );
result.at( i + half ) = Dif( lower1, upper1 );
o = Mul( o, oBase );
}
// My knowledge of move-semantics of STL containers is a bit rusty - so there's probably a better way to return the output 'result' vector.
return result;
}
FFT 工作正常,但当我想使用 IFFT 时,我总是从其结果中看到相同的图形。结果很复杂,无论原始信号如何,图形始终相同。
实部图是一个 -sin,周期 = 帧大小
虚部是同周期的-cos
哪里会出问题?
原始信号:
IFFT真实值(图片只有一半帧):
我用的算法FFT。
double** FFT(double** f, int s, bool inverse) {
if (s == 1) return f;
int sH = s / 2;
double** fOdd = new double*[sH];
double** fEven = new double*[sH];
for (int i = 0; i < sH; i++) {
int j = 2 * i;
fOdd[i] = f[j];
fEven[i] = f[j + 1];
}
double** sOdd = FFT(fOdd, sH, inverse);
double** sEven = FFT(fEven, sH, inverse);
double**spectr = new double*[s];
double arg = inverse ? DoublePI / s : -DoublePI / s;
double*oBase = new double[2]{ cos(arg),sin(arg) };
double*o = new double[2]{ 1,0 };
for (int i = 0; i < sH; i++) {
double* sO1 = Mul(o, sOdd[i]);
spectr[i] = Sum(sEven[i], sO1);
spectr[i + sH] = Dif(sEven[i], sO1);
o = Mul(o, oBase);
}
return spectr;
}
"butterfly" 部分错误地应用了系数:
for (int i = 0; i < sH; i++) {
double* sO1 = sOdd[i];
double* sE1 = Mul(o, sEven[i]);
spectr[i] = Sum(sO1, sE1);
spectr[i + sH] = Dif(sO1, sE1);
o = Mul(o, oBase);
}
旁注:
我保留了你的注释,但这让事情变得混乱:
fOdd
有索引 0, 2, 4, 6, ... 所以它应该是 fEven
fEven
有索引 1, 3, 5, 7, ... 所以它应该是 fOdd
实际上 sOdd
应该是 sLower
而 sEven
应该是 sUpper
因为它们对应于 0:s/2
和 s/2:s-1
元素光谱分别为:
sLower = FFT(fEven, sH, inverse); // fEven is 0, 2, 4, ...
sUpper = FFT(fOdd, sH, inverse); // fOdd is 1, 3, 5, ...
那么蝴蝶就变成了:
for (int i = 0; i < sH; i++) {
double* sL1 = sLower[i];
double* sU1 = Mul(o, sUpper[i]);
spectr[i] = Sum(sL1, sU1);
spectr[i + sH] = Dif(sL1, sU1);
o = Mul(o, oBase);
}
这样写比较容易与 pseudocode example on wikipedia.
进行比较@Dai 是正确的,你会泄漏大量内存
关于内存,您可以使用 std::vector
来封装动态分配的数组,并确保在执行离开范围时释放它们。您可以使用 unique_ptr<double[]>
但性能提升在 IMO 上不值得,并且您失去了 at()
方法的安全性。
(基于@Robb 的回答)
其他一些提示:
- 避免神秘的标识符 - 程序应该是可读的,像“
f
”和“s
”这样的名称会使您的程序更难阅读和维护。 - 基于类型的匈牙利符号不受欢迎,因为现代编辑器会自动显示类型信息,因此它给标识符名称增加了不必要的复杂性。
- 索引使用
size_t
,而不是int
- STL 是你的朋友,使用它!
- 通过使用
const
预防只读数据的意外突变来预防错误。
像这样:
#include <vector>
using namespace std;
vector<double> fastFourierTransform(const vector<double> signal, const bool inverse) {
if( signal.size() < 2 ) return signal;
const size_t half = signal.size() / 2;
vector<double> lower; lower.reserve( half );
vector<double> upper; upper.reserve( half );
bool isEven = true;
for( size_t i = 0; i < signal.size(); i++ ) {
if( isEven ) lower.push_back( signal.at( i ) );
else upper.push_back( signal.at( i ) );
isEven = !isEven;
}
vector<double> lowerFft = fastFourierTransform( lower, inverse );
vector<double> upperFft = fastFourierTransform( upper, inverse );
vector<double> result;
result.reserve( signal.size() );
double arg = ( inverse ? 1 : -1 ) * ( DoublePI / signal.size() );
// Ideally these should be local `double` values passed directly into `Mul`.
unique_ptr<double[]> oBase = make_unique<double[]>( 2 );
oBase[0] = cos(arg);
oBase[1] = sin(arg);
unique_ptr<double[]> o = make_unique<double[]>( 2 );
o[0] = 0;
o[1] = 0;
for( size_t i = 0; i < half; i++ ) {
double* lower1 = lower.at( i );
double* upper1 = Mul( o, upper.at( i ) );
result.at( i ) = Sum( lower1, upper1 );
result.at( i + half ) = Dif( lower1, upper1 );
o = Mul( o, oBase );
}
// My knowledge of move-semantics of STL containers is a bit rusty - so there's probably a better way to return the output 'result' vector.
return result;
}