24小时热门版块排行榜

返回列表

【奖励】本帖被评价1次，作者baobiao007增加金币 0.8 个

baobiao007

木虫 (职业作家)

应助: 201 (大学生)
金币: 6482.7
帖子: 3050
在线: 1009.9小时
虫号: 505962

[资源] 【分享】频率域镶边法带通滤波c程序【原创】

#include
#include
#include
#include"FFT.h"
//效果似乎不如时域褶积法
/*构建镶边带通滤波器
  f0--通频带中心频率
  hf--通频带半宽度
  num-滤波器长度
  df--频率采样间隔
*/
void Filter(float h[],int num,float f0,float hf,float df)
{
int i;
double pi=3.1415926;
float bf=hf;//镶边宽度
float f1,f2,f3,f4,f;
f2=f0-hf; f3=f0+hf;
f1=f2-bf; f4=f3+bf;
for(i=0; i<=num/2; i++)
{
f=i*df;
if(f else if(f>=f1 && f<=f2)
h[ i ]=pow(sin(pi*(f-f1)/(2.0*(f2-f1))),2);
else if(f>f2 && f h[ i ]=1.0;
else if(f>=f3 && f<=f4)
h[ i ]=pow(sin(pi*(f-f4)/(2.0*(f4-f3))),2);
else
h[ i ]=0.0;
}
for(i=num/2+1; i h[ i ]=h[num-i];
}

void main()
{
const int M=128;//滤波器长度
float dt=0.002;//时间采样间隔
float df=1.0/(M*dt);//频率采样间隔
float Inisigr[M]={0};//要滤波的信号
float Inisigi[M]={0};
float Afsig2r[M]={0};//滤波后的信号
float Afsig2i[M]={0};
float Afsig[M]={0};//理想情况下滤波后的信号
float hr[M],hi[M];//滤波器频谱
int i,k;
k=(int)(log(M)/log(2.0)+0.5);
FILE *fp1,*fp2,*fp3;
fp1=fopen("inisig.txt","w" );//合成信号
      fp2=fopen("afsig.txt","w" );//理想50hz信号
fp3=fopen("afsig2.txt","w" );//滤波得到的50hz信号

//产生信号,20hz,50hz,80hz正弦信号叠加
double pi=3.1415926;
float f1,f2,f3;
f1=20.0;
f2=50.0;
f3=80.0;
for(i=0; i { Afsig[ i ]=sin(2*pi*f2*dt*i);
Inisigr[ i ]=Afsig[ i ]+sin(2*pi*f1*i*dt)+sin(2*pi*i*f3*dt);
                           fprintf(fp1,"%f\n",Inisigr[ i ] );
}

//形成滤波器
Filter(hr,M,50.0,10.0,df);
//对要滤波的信号进行FFT
fft(Inisigr,Inisigi,k,1);
//频谱相乘,滤去低于40hz高于60hz的信号
for(i=0; i {
Afsig2r[ i ]=hr[ i ]*Inisigr[ i ];
Afsig2i[ i ]=Inisigi[ i ]*hr[ i ];
}
//FFT反变换
fft(Afsig2r,Afsig2i,k,-1);

//输出结果
      for(i=0; i {
                     fprintf(fp2,"%f\n",Afsig[ i ]);
                     fprintf(fp3,"%f\n",Afsig2r[ i ]);
}
            fclose(fp1);
fclose(fp2);
fclose(fp3);

}
滤波前：
滤波后：

[ Last edited by baobiao007 on 2011-2-14 at 22:55 ]

回复此楼