FFT filter

[zhiyayan]

我有一个很白痴的问题。就是我有一组数据，差不多就是X[1000000],Y[1000000]这样的数据，X是时间列，存的就是时间，0.001，0.002……Y里边是某一个物理量的数据。我想对数据做相当于origin里边的FFT filter，滤波到某一个频率，比如扔掉5Hz以上的。怎么做呢？求指教～最好给一段代码啦。。

[zhiyayan]

[Quote=引用 1 楼 greatliudy123 的回复:]

没做过，FFT，路过帮顶
[/Quote]

谢谢~

[lp_grace]

把物理量的数据做fft变换，小于5hz的去掉就行了，然后再做反变换，跟时间没有关系

[lp_grace]

这个网上到处都是啊

void FFT(complex<double> * TD, complex<double> * FD, int power)
{
int count;
int i,j,k,bfsize,p;
double angle;
complex<double> *W,*X1,*X2,*X;

/*计算傅立叶变换点数*/
count=1<<power;

/*分配运算所需存储器*/
W=(complex<double> *)malloc(sizeof(complex<double>)*count/2);
X1=(complex<double> *)malloc(sizeof(complex<double>)*count);
X2=(complex<double> *)malloc(sizeof(complex<double>)*count);

/*计算加权系数*/
for(i=0;i<count/2;i++)
{
angle=-i*pi*2/count;
W[i]=complex<double>(cos(angle),sin(angle));
}

/*将时域点写入存储器*/
memcpy(X1,TD,sizeof(complex<double>)*count);

/*蝶形运算*/
for(k=0;k<power;k++)
{
for(j=0;j<1<<k;j++)
{
bfsize=1<<(power-k);
for(i=0;i<bfsize/2;i++)
{
p=j*bfsize;
X2[i+p]=X1[i+p]+X1[i+p+bfsize/2];
X2[i+p+bfsize/2]=(X1[i+p]-X1[i+p+bfsize/2])*W[i*(1<<k)];
}
}
X=X1;
X1=X2;
X2=X;
}

/*重新排序*/
for(j=0;j<count;j++)
{
p=0;
for(i=0;i<power;i++)
{
if (j&(1<<i))
p+=1<<(power-i-1);
}
FD[j]=X1[p];
}

/*释放存储器*/
free(W);
free(X1);
free(X2);
}

void IFFT(complex<double> * FD, complex<double> * TD, int power)
{
int i;
long count;
complex<double> *X;

/*计算傅立叶反变换点数*/
count=1<<power;

/*分配运算所需存储器*/
X=(complex<double> *)malloc(sizeof(complex<double>)*count);

/*将频域点写入存储器*/
memcpy(X,FD,sizeof(complex<double>)*count);

/*求频域点的共轭*/
for(i=0;i<count;i++)
X[i]=complex<double>(X[i].real(),-X[i].imag());

/*调用FFT*/
FFT(X, TD, power);

/*求时域点的共轭*/
for(i=0;i<count;i++)
TD[i]=complex<double>(TD[i].real()/count,-TD[i].imag()/count);

/*释放存储器*/
free(X);
}[Quote=引用 6 楼 zhiyayan 的回复:]

引用 5 楼 lp_grace 的回复:

fft(X,Y);//X时域数据，Y频域数据
for（int i=0;i<n;i++）
{
if(Y[i]<5)
Y[i]=0;
}

引用 4 楼 zhiyayan 的回复:

引用 3 楼 lp_grace 的回复:

把物理量的数据做fft变换，小于5hz的去掉就行了，然后再做反变换，跟时间没有关系
……
[/Quote]

[zhiyayan]

非常感谢你！然后TD,FD,power分别是什么呢？power怎么设置这个参数呢？还有我的滤波频率f怎么和FFT变换之后数组中的数对应呢？比如滤波到10.0223Hz，第几个数才是相对的滤波频率呢？也就是说从那个数开始置零？真不好意思，这就是看算法所不明白的地方。。[Quote=引用 7 楼 lp_grace 的回复:]

这个网上到处都是啊

void FFT(complex<double> * TD, complex<double> * FD, int power)
{
int count;
int i,j,k,bfsize,p;
double angle;
complex<double> *W,*X1,*X2,*X;

/*计算傅立叶变换点数*/
coun……
[/Quote]

[zhiyayan]

如果没有时间轴列，怎么定义频率呢？
[Quote=引用 7 楼 lp_grace 的回复:]

这个网上到处都是啊

void FFT(complex<double> * TD, complex<double> * FD, int power)
{
int count;
int i,j,k,bfsize,p;
double angle;
complex<double> *W,*X1,*X2,*X;

/*计算傅立叶变换点数*/
coun……
[/Quote]

[zhiyayan]

[Quote=引用 3 楼 lp_grace 的回复:]

把物理量的数据做fft变换，小于5hz的去掉就行了，然后再做反变换，跟时间没有关系
[/Quote]

这个我知道是去掉5HZ以上的。然后关键是怎么做？看到教程里的算法根本不明白是怎么一回事。。求代码。。

[lp_grace]

fft(X,Y);//X时域数据，Y频域数据
for（int i=0;i<n;i++）
{
if(Y[i]<5)
Y[i]=0;
}

[Quote=引用 4 楼 zhiyayan 的回复:]

引用 3 楼 lp_grace 的回复:

把物理量的数据做fft变换，小于5hz的去掉就行了，然后再做反变换，跟时间没有关系


这个我知道是去掉5HZ以上的。然后关键是怎么做？看到教程里的算法根本不明白是怎么一回事。。求代码。。
[/Quote]

[zhiyayan]

[Quote=引用 5 楼 lp_grace 的回复:]

fft(X,Y);//X时域数据，Y频域数据
for（int i=0;i<n;i++）
{
if(Y[i]<5)
Y[i]=0;
}

引用 4 楼 zhiyayan 的回复:

引用 3 楼 lp_grace 的回复:

把物理量的数据做fft变换，小于5hz的去掉就行了，然后再做反变换，跟时间没有关系


这个我知道是去掉5HZ以上的。然后关键是……
[/Quote]

我说的是FFT和IFFT的代码呃。。谢谢啦～