24小时热门版块排行榜    

查看: 1018  |  回复: 1
【奖励】 本帖被评价1次,作者baobiao007增加金币 0.8

baobiao007

木虫 (职业作家)


[资源] 【分享】频率域镶边法带通滤波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 ]
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

w1161392713

新虫 (初入文坛)


★★★★★ 五星级,优秀推荐

请问有没有这个快速傅里叶变换的头文件FFT.h?
2楼2017-05-20 13:45:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 baobiao007 的主题更新
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 投票:  有多少人是今天查系统知道结果的? +17 爱看书的可乐 2026-08-26 19/950 2026-09-02 00:19 by xiangy672
[基金申请] 面上函评意见出来了,像什么等级? 20+4 Tsingking1 2026-08-27 17/850 2026-09-01 19:51 by 超级无敌华子
[基金申请] 为什么国自然不能直接公布 +5 bjdxyxy 2026-08-26 5/250 2026-09-01 17:54 by 小伟大博士
[文学芳草园] 梦想 +5 myrtle 2026-08-26 7/350 2026-09-01 15:18 by myrtle
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +7 gdd2018 2026-08-28 12/600 2026-09-01 08:32 by 尼古拉斯小虫
[基金申请] 基金未中,这种答复是模板吗? +6 zhaosm1982 2026-08-27 7/350 2026-08-31 21:18 by qdxxmc
[基金申请] 面上意见出来了 +12 黄鸟于飞Chao 2026-08-29 23/1150 2026-08-31 18:57 by 黄鸟于飞Chao
[基金申请] 能否申诉? +7 echo8914667 2026-08-30 8/400 2026-08-31 17:00 by yihongxu
[基金申请] 中青基了要发朋友圈吗? +7 349506619 2026-08-28 7/350 2026-08-31 13:39 by 冼亮淀粉酶
[基金申请] 有没有仍没收到信息的 +7 德尚中行 2026-08-27 8/400 2026-08-30 20:52 by purplejack
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 国自然面上复盘~欢迎讨论 (金币+15) +15 晴天加油 2026-08-26 16/800 2026-08-29 18:28 by symmetry
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 我不理解! +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 国合里面能看到了 +7 一怀馨秋 2026-08-26 7/350 2026-08-26 11:23 by zhaosm1982
[基金申请] 国际合作可查了,中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 国合可查了 +3 paperzjh 2026-08-26 3/150 2026-08-26 10:41 by LemmonTr
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
信息提示
请填处理意见