24小时热门版块排行榜    

查看: 1014  |  回复: 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 的主题更新
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] 学科评审组评审是指会评吗? +4 瞬息宇宙 2026-08-31 4/200 2026-09-01 14:58 by jiaoxg
[基金申请] 面上合作单位盖章 +5 ssyjh 2026-08-27 8/400 2026-09-01 12:58 by ssyjh
[基金申请] 面上函评意见出来了,像什么等级? 20+4 Tsingking1 2026-08-27 15/750 2026-09-01 12:23 by icm639
[基金申请] 麻烦专家们看看评委们的意见(F口面上) +7 gdd2018 2026-08-28 12/600 2026-09-01 08:32 by 尼古拉斯小虫
[基金申请] 怎么看青基中了没有啊 +6 叶九微 2026-08-26 6/300 2026-08-31 23:54 by yudaoqian88
[基金申请] 国社科又开始会评了,不知道这次命运如何 +7 雨打竹帘 2026-08-30 11/550 2026-08-31 23:16 by hittle2008
[基金申请] 投票:  有多少人是今天查系统知道结果的? +17 爱看书的可乐 2026-08-26 19/950 2026-08-31 21:30 by xiangy672
[基金申请] 基金未中,这种答复是模板吗? +6 zhaosm1982 2026-08-27 7/350 2026-08-31 21:18 by qdxxmc
[基金申请] 基金不中,共勉 +12 eulota 2026-08-26 12/600 2026-08-31 08:42 by ZJTJZ
[基金申请] 为什么到现在没收到通知? +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[基金申请] 有没有仍没收到信息的 +7 德尚中行 2026-08-27 8/400 2026-08-30 20:52 by purplejack
[考博] 找导师 +6 yuanjiabao 2026-08-29 7/350 2026-08-30 14:40 by 生科新手
[基金申请] 2026年叶企孙基金 +4 bud_bud 2026-08-27 7/350 2026-08-29 07:23 by foolishmani
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[基金申请] 我不理解! +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
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 项目信息和经费信息在系统里都可以看到了 +6 wittyboy 2026-08-26 14/700 2026-08-26 10:55 by wittyboy
[基金申请] 牛来!米来!面来! +8 beefly 2026-08-26 8/400 2026-08-26 08:37 by xuzhipiao
信息提示
请填处理意见