文章
目录
前言
现在大多教程都是matlab的教程，那我用C语言实现一下FIR
此博客为个人学习
一、FIR滤波器
FIR有限长滤波器，有限长指的是滤波器的脉冲响应为有限长的序列
输入X进过滤波器H得到滤波后的输出Y
表达成差分方程形式就是：
也可以写成卷鸡形式
y[n]=h[n]*x[n]
窗函数法
由于理想滤波器在时域上是无线长序列，所以要对理想滤波器加窗截断变成有限信号。
理想低通滤波器时域
常用窗函数：
使用不同的窗函数会有不同的滤波性能（衰减程度和过渡带大小的不一样）
窗函数的介绍以及画出常见窗函数（汉宁窗，矩形窗，汉明窗，布莱克曼窗）的时域图和频谱图-CSDN博客
设计步骤
构造理想滤波器
hl[n]
理想离散低通滤波器公式
低通可构成其他高通带通等
构造窗函数
w[n]
相乘得到FIR滤波器
h[n]=w[n]hl[n]
最后让输出
得到滤波数据
二、代码
0.头文件和宏定义
#include<stdio.h>
#include <stdlib.h>
#include <math.h>
#define pi 3.1415926535897932
/*----窗类型-----*/
//三角形窗未完成，所以其实三角形和矩形的一样
#define Rectangle 1 //矩形
#define triangle 2 //三角
#define Hanning 3 //汉宁
#define Hamming 4 //海明
#define Blackman 5 //布拉克曼
/*----滤波类型-----*/
#define LOWPASSFILTER 1 //低通
#define HIGEPASSFILTER 2
#define BANDPASSFILTER 3 //带通
#define BANDSTOPFILTER 41.构造滤波器1函数
这个函数是按阶数设计滤波器,(N为偶数时，窗函数法无法设计高通和带组)
//已知阶数，求fir 的 传递函数 值
/*-----------------------
* 函数名 FIR_Filter_Transfer_functions_param
* 输入 @Filter_h ： 用于返回传递函数 数据数组 长度为N
@Filter_type ：滤波器类型 低通高通带通带阻 参数可以是[LOWPASSFILTER , HIGEPASSFILTER ,BANDPASSFILTER , BANDSTOPFILTER]
@Wc1: 滤波器的截止频率
@Wc2： BPF和BSF的第二个截止频率
@N 阶数
@window 窗函数类型：参数可以是[Rectangle, triangle , Hanning ,Hamming ,Blackman]
* 功能 按阶数设计滤波器 ，注意设计高通和带阻时 N必须为奇数
* 返回 1成功 0失败
----------------------*/
char FIR_Filter_Transfer_functions_param(double *Filter_h,char Filter_type,double Wc1 ,double Wc2 ,unsigned int N,int window)
{
int n;
//偶数点无法设计高通和带阻 //
if(((Filter_type==HIGEPASSFILTER)||(Filter_type==BANDSTOPFILTER))&&(N%2==0))
{
return 0;
}
/*窗函数******-------------------------------*/
int const const_N=N;
double win_param[const_N]; //存放窗函数数据
switch(window)
{
case Rectangle:
for(n=0;n<N;n++)
{
win_param[n]=1;
}
break;
case triangle:
for(n=0;n<N;n++)
{
win_param[n]=1; //2部分，懒得打了 ，直接用矩形窗游览
}
break;
case Hanning:
for(n=0;n<N;n++)
{
win_param[n]=0.5*(1-cos(2*pi*n/(N-1)));
}
break;
case Hamming :
for(n=0;n<N;n++)
{
win_param[n]=0.54-0.46*cos(2*pi*n/(N-1));
// printf("%lf\n",win_param[n]);
}
break;
case Blackman :
for(n=0;n<N;n++)
{
win_param[n]=0.42-0.5*cos((4*pi*n/(N-1))-0.08*(4*pi*n/(N-1)));
}
break;
default:
return 0;
}
/*---------------构造理想滤波器-------------------------*/
double tao; //群延迟
double h_lixiang[const_N]; //理想滤波器
tao=(N-1)/2.0;
switch(Filter_type)
{
case LOWPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //如果 (n-tao)==0则分子分母为0 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=sin(Wc1*(n-tao))/(pi*(n-tao));
}
}
break;
case HIGEPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0)
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(pi*(n-tao))-sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
case BANDPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0)
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(Wc2*(n-tao))-sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
case BANDSTOPFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0)
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(pi*(n-tao))-sin(Wc2*(n-tao))+sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
default:
return 0;
}
/*--------求h[n]-----------*/
for(n=0;n<N;n++)
{
Filter_h[n]=h_lixiang[n]*win_param[n];
}
return 1;
}
2.构造滤波器2函数
这个是按照模拟指标设计滤波器
配套结构体
typedef struct{
double(*Hn)[]; // 指向一个 传递函数的指针
int N; //传递函数的阶数
}H_Struct; /*传递函数结构体------第二种函数实现用到 */
typedef struct {
int fs; //采样率
int fp1; //通带截止频率 1 靠近w=0这边的
int fst1; //通带截止频率 1 靠近w=0这边的
int fp2; //通带截止频率 2 如果为低通和高通则 不用关 2
int fst2; //阻带截止频率 2
int ast; //阻带最小衰减
}Transfer_functions_struct; /*模拟指标结构体 */
/*-----------------------
* 函数名 FIR_Filter_Transfer_functions_param2
* 输入
@Filter_type ： 滤波器类型 低通高通带通带阻 参数可以是[LOWPASSFILTER , HIGEPASSFILTER ,BANDPASSFILTER , BANDSTOPFILTER]
@Transfer_param: 滤波器的模拟指标，
* 功能 按照模拟指标设计滤波器
* 返回 H_Struct包含 滤波器传递函数 和 阶数 的结构体 ，
----------------------*/
H_Struct* FIR_Filter_Transfer_functions_param2(char Filter_type,Transfer_functions_struct*Transfer_param)
{
/*------------数字指标-------------------*/
double Wc1=2*pi*((Transfer_param->fp1+Transfer_param->fst1)/2)/Transfer_param->fs;
double Wc2=2*pi*((Transfer_param->fp2+Transfer_param->fst2)/2)/Transfer_param->fs;
double det_W1=2*pi*abs(Transfer_param->fst1-Transfer_param->fp1)/Transfer_param->fs; //过度带带宽
double det_W2=2*pi*abs(Transfer_param->fst2-Transfer_param->fp2)/Transfer_param->fs; //过度带带宽 2
double det_W=det_W1;
if((Filter_type==BANDPASSFILTER)||(Filter_type==BANDSTOPFILTER)) //如果是带通或者带阻，选择窄的指标
{
det_W=det_W1>det_W2?det_W2:det_W1;
}
int n;
int N;
int window;
/*窗函数选择--------------------------------------------------------------------------- */
// 根据衰as减选择窗
switch((int)(Transfer_param->ast/10))
{
case 2: //20db
N=(int)1.8*pi/det_W;
if(N!=1.8*pi/det_W)N++;
window=Rectangle;
break;
case 3://30db
N=(int)6.1*pi/det_W;
if(N!=6.1*pi/det_W)N++;
window=triangle;
break;
case 4: //40db
N=(int)(6.2*pi/det_W);
if(N!=6.2*pi/det_W)N++;
window=Hanning;
break;
case 5 ://50db
N=(int)(6.6*pi/det_W);
if(N!=6.6*pi/det_W)N++;
window=Hamming;
break;
case 7 ://70db
N=(int)11*pi/det_W;
if(N!=11*pi/det_W)N++;
window=Blackman;
break;
default:
return 0;
}
int const const_N=N;
double win_param[const_N]; //存放窗函数数据
switch(window)
{
case Rectangle: //矩形
for(n=0;n<N;n++)
{
win_param[n]=1;
}
break;
case triangle:
for(n=0;n<N;n++)
{
win_param[n]=1; //三角形窗未实现
}
break;
case Hanning:
for(n=0;n<N;n++)
{
win_param[n]=0.5*(1-cos(2*pi*n/(N-1)));
}
break;
case Hamming :
for(n=0;n<N;n++)
{
win_param[n]=0.54-0.46*cos(2*pi*n/(N-1));
// printf("%lf\n",win_param[n]);
}
break;
case Blackman :
for(n=0;n<N;n++)
{
win_param[n]=0.42-0.5*cos((4*pi*n/(N-1))-0.08*(4*pi*n/(N-1)));
}
break;
default: //只设计了几个常用的，衰减超出这些则返回失败
return 0;
}
/*构造理想滤波器-----------------------------------------------------------------------------*/
double tao; //H的群延迟
double h_lixiang[const_N]; //存放理想滤波器 数据
tao=(N-1)/2.0; //求出中心点
switch(Filter_type)
{
case LOWPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //n-tao=0时分子分母为零 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=sin(Wc1*(n-tao))/(pi*(n-tao)); //理想低通滤波器表达式
}
}
break;
case HIGEPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //n-tao=0时分子分母为零 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(pi*(n-tao))-sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
case BANDPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //n-tao=0时分子分母为零 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(Wc2*(n-tao))-sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
case BANDSTOPFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //n-tao=0时分子分母为零 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(pi*(n-tao))-sin(Wc2*(n-tao))+sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
default: //如果不是这些则返回失败
return 0;
}
/*计算FIR的传递函数*/
H_Struct *Filter_h=NULL;
Filter_h=(H_Struct*)malloc(sizeof(H_Struct));
Filter_h->Hn=(double(*)[])malloc(sizeof(double)*N); //申请 跟阶数相等的内存空间
Filter_h->N=N;
for(n=0;n<N;n++) //将窗函数和理想低通相乘得到 h
{
(*(Filter_h->Hn))[n]=h_lixiang[n]*win_param[n];
}
return Filter_h;
}3滤波函数
配套结构体和初始化函数
typedef struct {
double(*input_Xbuff)[]; //指向一个缓存区用于存放历史数据 ，长度为N的数组
int Data_input;
int Data_output ;
}FIR_struct; /*一个存历史输入和当前输入的 环形数组 */
*-----------------------
* 函数名 FIR_Struct_Init
* 输入 @Fir_variable ：一个 历史输入 环形数组指针
* 功能 初始化数组
* 返回 无
----------------------*/
void FIR_Struct_Init(FIR_struct *Fir_variable)
{
int i;
Fir_variable->Data_input=0;
Fir_variable->Data_output=0;
Fir_variable->input_Xbuff=NULL;
}
函数
/*-----------------------
* 函数名 Fir_filter
* 输入 @Xdata： 包含这一时刻的 历史数据数组
@hn ： 滤波器的传递函数
@N 阶数
* 功能 求出滤波后的输出
* 返回 这一时刻滤波后的输出
----------------------*/
double Fir_filter(FIR_struct*Xdata,double*hn,unsigned int N)
{
double y=0;
int i=0,j=0;
for(i=0;i<N;i++)
{
//寻找历史输入的数组 下标 X[n-i]
if((Xdata->Data_output-i)<0)
{
//越界
j=N+Xdata->Data_output-i;
}else
{
j=Xdata->Data_output-i;
}
/*Fir的差分方程形式-->卷积*/
y+=hn[i]*(*(Xdata->input_Xbuff))[j];
}
/*延时数据,此刻变成上一时刻*/
Xdata->Data_output++;
if(Xdata->Data_output>=N)
{
Xdata->Data_output=0;
}
return y;
}三、使用步骤和matlab验证
输入信号是320 、1200 、2400hz叠加的波形
double X[500];
double y;
/*---*/
double old[Nd]; //历史输入数组
FIR_struct X_old;
FIR_Struct_Init(&X_old);
X_old.input_Xbuff=&old;
int t=500;
/*---输入信号*/
for(k=0;k<t;k++)
{
X[k]=sin(2*pi*320*k/fs)+sin(2*pi*1200*k/fs)+sin(2*pi*2400*k/fs);
}
1.用阶数设计滤波器
设计一个通带Nd=42阶 的低通滤波 截止频率fc在1000 到1400hz之间的低通滤波
double hn[Nd]={0}; //fir滤波器数组
double fs=5000;
double fp=1000;
double fst=1400;
double wc=2*pi*((fst+fp)/2)/fs; //3db频率
int i=0,k;
/*---*/
double old[Nd]; //历史输入数组
FIR_struct X_old;
FIR_Struct_Init(&X_old);
X_old.input_Xbuff=&old;
FIR_Filter_Transfer_functions_param(hn,LOWPASSFILTER,wc,0,Nd,Hamming);
FILE*H_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/H_pm.txt", "w"); //输出H的数据文件让matlab调用绘图
for(i=0;i<Nd;i++)
{
fprintf(H_fp, "%lf\n",hn[i]);
}
FILE*Y_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/Y_pm.txt", "w"); //输出Y的数据文件让matlab调用绘图
for(i=0;i<t;i++)
{
(*(X_old.input_Xbuff))[X_old.Data_input]=X[i];
X_old.Data_input++;
X_old.Data_input%=Nd;
y=Fir_filter(&X_old,hn,Nd);
fprintf(Y_fp, "%lf\n",y); //写入数据到输出文件数据
}
fclose(Y_fp);//关闭文件
matlab绘制数据图形，可以看出滤波效果不错
（滤波前后对比）
（滤波器幅频响应）
2.用指标设计滤波器
设计一个带通滤波器
指标：
采样频率=5000
指标as=50
通带截止为800 和1500
阻带截止为400 和2000
Transfer_functions_struct H;
H.ast=50;
H.fp1=800;
H.fst1=400;
H.fp2=1500;
H.fst2=2000;
H.fs=5000;
H_Struct *hk=NULL;
hk=FIR_Filter_Transfer_functions_param2(BANDSTOPFILTER,&H);
FILE*H_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/H_pm.txt", "w"); //输出H的数据文件让matlab调用
for(i=0;i<hk->N;i++)
{
fprintf(H_fp, "%lf\n",(*(hk->Hn))[i]);
}
fclose(H_fp);
FILE*Y_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/Y_pm.txt", "w"); //输出Y的数据文件让matlab调用
for(i=0;i<t;i++)
{
(*(X_old.input_Xbuff))[X_old.Data_input]=X[i]; //将 输入放入历史输入数组
X_old.Data_input++;
X_old.Data_input%=hk->N;
y=Fir_filter(&X_old,(*(hk->Hn)),hk->N);
//printf("%lf\n",y);
fprintf(Y_fp, "%lf\n",y);
}
fclose(Y_fp);
printf("%d\n",hk->N);matlab
（滤波器幅频响应）
四、个人实验的完整代码
#include<stdio.h>
#include <stdlib.h>
#include <math.h>
#define pi 3.1415926535897932
/*----窗类型-----*/
//三角形窗未完成，所以其实三角形和矩形的一样
#define Rectangle 1 //矩形
#define triangle 2 //三角
#define Hanning 3 //汉宁
#define Hamming 4 //海明
#define Blackman 5 //布拉克曼
/*----滤波类型-----*/
#define LOWPASSFILTER 1 //低通
#define HIGEPASSFILTER 2
#define BANDPASSFILTER 3 //带通
#define BANDSTOPFILTER 4
typedef struct {
double(*input_Xbuff)[]; //指向一个缓存区用于存放历史数据 ，长度为N的数组
int Data_input;
int Data_output ;
}FIR_struct; /*一个存历史输入和当前输入的 环形数组 */
typedef struct{
double(*Hn)[]; // 指向一个 传递函数的指针
int N; //传递函数的阶数
}H_Struct; /*传递函数结构体------第二种函数实现用到 */
typedef struct {
int fs; //采样率
int fp1; //通带截止频率 1 靠近w=0这边的
int fst1; //通带截止频率 1 靠近w=0这边的
int fp2; //通带截止频率 2 如果为低通和高通则 不用关 2
int fst2; //阻带截止频率 2
int ast; //阻带最小衰减
}Transfer_functions_struct; /*模拟指标结构体 */
void FIR_Struct_Init(FIR_struct *Fir_variable);
H_Struct* FIR_Filter_Transfer_functions_param2(char Filter_type,Transfer_functions_struct*Transfer_param);
char FIR_Filter_Transfer_functions_param(double *Filter_h,char Filter_type,double Wc1 ,double Wc2 ,unsigned int N,int window);
double Fir_filter(FIR_struct*Xdata,double*hn,unsigned int N);
#define Nd 42 //阶数
#define Typ 3 //当 1 2 不同函数设计低通 否则用第二种函数设计 带通
int main()
{
double hn[Nd]={0}; //fir滤波器数组
double fs=5000;
double fp=1000;
double fst=1400;
double wc=2*pi*((fst+fp)/2)/fs; //3db频率
int i=0,k;
double X[500];
double y;
/*---*/
double old[Nd]; //历史输入数组
FIR_struct X_old;
FIR_Struct_Init(&X_old);
X_old.input_Xbuff=&old;
int t=500;
/*---输入信号*/
for(k=0;k<t;k++)
{
X[k]=sin(2*pi*320*k/fs)+sin(2*pi*1200*k/fs)+sin(2*pi*2400*k/fs);
}
#if Typ==1 /*用第一种函数设计低通滤波器*/
FIR_Filter_Transfer_functions_param(hn,LOWPASSFILTER,wc,0,Nd,Hamming);
FILE*H_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/H_pm.txt", "w"); //输出H的数据文件让matlab调用绘图
for(i=0;i<Nd;i++)
{
fprintf(H_fp, "%lf\n",hn[i]);
}
FILE*Y_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/Y_pm.txt", "w"); //输出Y的数据文件让matlab调用绘图
for(i=0;i<t;i++)
{
(*(X_old.input_Xbuff))[X_old.Data_input]=X[i];
X_old.Data_input++;
X_old.Data_input%=Nd;
y=Fir_filter(&X_old,hn,Nd);
fprintf(Y_fp, "%lf\n",y); //写入数据到输出文件数据
}
fclose(Y_fp);//关闭文件
#elif Typ==2 /*用第二种函数设计低通滤波器*/
Transfer_functions_struct H;
H.ast=50;
H.fp1=1000;
H.fst1=1400;
H.fs=5000;
H_Struct *hk=NULL;
hk=FIR_Filter_Transfer_functions_param2(LOWPASSFILTER,&H);
FILE*H_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/H_pm.txt", "w"); //输出H的数据文件让matlab调用
for(i=0;i<hk->N;i++)
{
fprintf(H_fp, "%lf\n",(*(hk->Hn))[i]);
}
fclose(H_fp);
FILE*Y_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/Y_pm.txt", "w"); //输出Y的数据文件让matlab调用
for(i=0;i<t;i++)
{
(*(X_old.input_Xbuff))[X_old.Data_input]=X[i];
X_old.Data_input++;
X_old.Data_input%=hk->N;
y=Fir_filter(&X_old,(*(hk->Hn)),hk->N);
fprintf(Y_fp, "%lf\n",y);
}
fclose(Y_fp);
printf("%d\n",hk->N);
#else /*设计一个带通滤波器*/
Transfer_functions_struct H;
H.ast=50;
H.fp1=800;
H.fst1=400;
H.fp2=1500;
H.fst2=2000;
H.fs=5000;
H_Struct *hk=NULL;
hk=FIR_Filter_Transfer_functions_param2(BANDSTOPFILTER,&H);
FILE*H_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/H_pm.txt", "w"); //输出H的数据文件让matlab调用
for(i=0;i<hk->N;i++)
{
fprintf(H_fp, "%lf\n",(*(hk->Hn))[i]);
}
fclose(H_fp);
FILE*Y_fp = fopen("E:/Matlab1/bin/DSP_poj/Fir_filter/Y_pm.txt", "w"); //输出Y的数据文件让matlab调用
for(i=0;i<t;i++)
{
(*(X_old.input_Xbuff))[X_old.Data_input]=X[i]; //将 输入放入历史输入数组
X_old.Data_input++;
X_old.Data_input%=hk->N;
y=Fir_filter(&X_old,(*(hk->Hn)),hk->N);
//printf("%lf\n",y);
fprintf(Y_fp, "%lf\n",y);
}
fclose(Y_fp);
printf("%d\n",hk->N);
#endif
return 0;
}
/*-----------------------
* 函数名 FIR_Struct_Init
* 输入 @Fir_variable ：一个 历史输入 环形数组指针
* 功能 初始化数组
* 返回 无
----------------------*/
void FIR_Struct_Init(FIR_struct *Fir_variable)
{
int i;
Fir_variable->Data_input=0;
Fir_variable->Data_output=0;
Fir_variable->input_Xbuff=NULL;
}
/*-----------------------
* 函数名 FIR_Filter_Transfer_functions_param2
* 输入
@Filter_type ： 滤波器类型 低通高通带通带阻 参数可以是[LOWPASSFILTER , HIGEPASSFILTER ,BANDPASSFILTER , BANDSTOPFILTER]
@Transfer_param: 滤波器的模拟指标，
* 功能 按照模拟指标设计滤波器
* 返回 H_Struct包含 滤波器传递函数 和 阶数 的结构体 ，
----------------------*/
H_Struct* FIR_Filter_Transfer_functions_param2(char Filter_type,Transfer_functions_struct*Transfer_param)
{
/*------------数字指标-------------------*/
double Wc1=2*pi*((Transfer_param->fp1+Transfer_param->fst1)/2)/Transfer_param->fs;
double Wc2=2*pi*((Transfer_param->fp2+Transfer_param->fst2)/2)/Transfer_param->fs;
double det_W1=2*pi*abs(Transfer_param->fst1-Transfer_param->fp1)/Transfer_param->fs; //过度带带宽
double det_W2=2*pi*abs(Transfer_param->fst2-Transfer_param->fp2)/Transfer_param->fs; //过度带带宽 2
double det_W=det_W1;
if((Filter_type==BANDPASSFILTER)||(Filter_type==BANDSTOPFILTER)) //如果是带通或者带阻，选择窄的指标
{
det_W=det_W1>det_W2?det_W2:det_W1;
}
int n;
int N;
int window;
/*窗函数选择--------------------------------------------------------------------------- */
// 根据衰as减选择窗
switch((int)(Transfer_param->ast/10))
{
case 2: //20db
N=(int)1.8*pi/det_W;
if(N!=1.8*pi/det_W)N++;
window=Rectangle;
break;
case 3://30db
N=(int)6.1*pi/det_W;
if(N!=6.1*pi/det_W)N++;
window=triangle;
break;
case 4: //40db
N=(int)(6.2*pi/det_W);
if(N!=6.2*pi/det_W)N++;
window=Hanning;
break;
case 5 ://50db
N=(int)(6.6*pi/det_W);
if(N!=6.6*pi/det_W)N++;
window=Hamming;
break;
case 7 ://70db
N=(int)11*pi/det_W;
if(N!=11*pi/det_W)N++;
window=Blackman;
break;
default:
return 0;
}
int const const_N=N;
double win_param[const_N]; //存放窗函数数据
switch(window)
{
case Rectangle: //矩形
for(n=0;n<N;n++)
{
win_param[n]=1;
}
break;
case triangle:
for(n=0;n<N;n++)
{
win_param[n]=1; //三角形窗未实现
}
break;
case Hanning:
for(n=0;n<N;n++)
{
win_param[n]=0.5*(1-cos(2*pi*n/(N-1)));
}
break;
case Hamming :
for(n=0;n<N;n++)
{
win_param[n]=0.54-0.46*cos(2*pi*n/(N-1));
// printf("%lf\n",win_param[n]);
}
break;
case Blackman :
for(n=0;n<N;n++)
{
win_param[n]=0.42-0.5*cos((4*pi*n/(N-1))-0.08*(4*pi*n/(N-1)));
}
break;
default: //只设计了几个常用的，衰减超出这些则返回失败
return 0;
}
/*构造理想滤波器-----------------------------------------------------------------------------*/
double tao; //H的群延迟
double h_lixiang[const_N]; //存放理想滤波器 数据
tao=(N-1)/2.0; //求出中心点
switch(Filter_type)
{
case LOWPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //n-tao=0时分子分母为零 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=sin(Wc1*(n-tao))/(pi*(n-tao)); //理想低通滤波器表达式
}
}
break;
case HIGEPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //n-tao=0时分子分母为零 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(pi*(n-tao))-sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
case BANDPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //n-tao=0时分子分母为零 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(Wc2*(n-tao))-sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
case BANDSTOPFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //n-tao=0时分子分母为零 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(pi*(n-tao))-sin(Wc2*(n-tao))+sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
default: //如果不是这些则返回失败
return 0;
}
/*计算FIR的传递函数*/
H_Struct *Filter_h=NULL;
Filter_h=(H_Struct*)malloc(sizeof(H_Struct));
Filter_h->Hn=(double(*)[])malloc(sizeof(double)*N); //申请 跟阶数相等的内存空间
Filter_h->N=N;
for(n=0;n<N;n++) //将窗函数和理想低通相乘得到 h
{
(*(Filter_h->Hn))[n]=h_lixiang[n]*win_param[n];
}
return Filter_h;
}
/*-----------------------
* 函数名 Fir_filter
* 输入 @Xdata： 包含这一时刻的 历史数据数组
@hn ： 滤波器的传递函数
@N 阶数
* 功能 求出滤波后的输出
* 返回 这一时刻滤波后的输出
----------------------*/
double Fir_filter(FIR_struct*Xdata,double*hn,unsigned int N)
{
double y=0;
int i=0,j=0;
for(i=0;i<N;i++)
{
//寻找历史输入的数组 下标 X[n-i]
if((Xdata->Data_output-i)<0)
{
//越界
j=N+Xdata->Data_output-i;
}else
{
j=Xdata->Data_output-i;
}
/*Fir的差分方程形式-->卷积*/
y+=hn[i]*(*(Xdata->input_Xbuff))[j];
}
/*延时数据,此刻变成上一时刻*/
Xdata->Data_output++;
if(Xdata->Data_output>=N)
{
Xdata->Data_output=0;
}
return y;
}
//已知阶数，求fir 的 传递函数 值
/*-----------------------
* 函数名 FIR_Filter_Transfer_functions_param
* 输入 @Filter_h ： 用于返回传递函数 数据数组 长度为N
@Filter_type ：滤波器类型 低通高通带通带阻 参数可以是[LOWPASSFILTER , HIGEPASSFILTER ,BANDPASSFILTER , BANDSTOPFILTER]
@Wc1: 滤波器的截止频率
@Wc2： BPF和BSF的第二个截止频率
@N 阶数
@window 窗函数类型：参数可以是[Rectangle, triangle , Hanning ,Hamming ,Blackman]
* 功能 按阶数设计滤波器 ，注意设计高通和带阻时 N必须为奇数
* 返回 1成功 0失败
----------------------*/
char FIR_Filter_Transfer_functions_param(double *Filter_h,char Filter_type,double Wc1 ,double Wc2 ,unsigned int N,int window)
{
int n;
//偶数点无法设计高通和带阻 //
if(((Filter_type==HIGEPASSFILTER)||(Filter_type==BANDSTOPFILTER))&&(N%2==0))
{
return 0;
}
/*窗函数******-------------------------------*/
int const const_N=N;
double win_param[const_N]; //存放窗函数数据
switch(window)
{
case Rectangle:
for(n=0;n<N;n++)
{
win_param[n]=1;
}
break;
case triangle:
for(n=0;n<N;n++)
{
win_param[n]=1; //2部分，懒得打了 ，直接用矩形窗游览
}
break;
case Hanning:
for(n=0;n<N;n++)
{
win_param[n]=0.5*(1-cos(2*pi*n/(N-1)));
}
break;
case Hamming :
for(n=0;n<N;n++)
{
win_param[n]=0.54-0.46*cos(2*pi*n/(N-1));
// printf("%lf\n",win_param[n]);
}
break;
case Blackman :
for(n=0;n<N;n++)
{
win_param[n]=0.42-0.5*cos((4*pi*n/(N-1))-0.08*(4*pi*n/(N-1)));
}
break;
default:
return 0;
}
/*---------------构造理想滤波器-------------------------*/
double tao; //群延迟
double h_lixiang[const_N]; //理想滤波器
tao=(N-1)/2.0;
switch(Filter_type)
{
case LOWPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0) //如果 (n-tao)==0则分子分母为0 取极限
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=sin(Wc1*(n-tao))/(pi*(n-tao));
}
}
break;
case HIGEPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0)
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(pi*(n-tao))-sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
case BANDPASSFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0)
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(Wc2*(n-tao))-sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
case BANDSTOPFILTER:
for(n=0;n<N;n++)
{
if((n-tao)==0)
{
h_lixiang[n]=Wc1/pi;
}
else
{
h_lixiang[n]=(sin(pi*(n-tao))-sin(Wc2*(n-tao))+sin(Wc1*(n-tao)))/(pi*(n-tao));
}
}
break;
default:
return 0;
}
/*--------求h[n]-----------*/
for(n=0;n<N;n++)
{
Filter_h[n]=h_lixiang[n]*win_param[n];
}
return 1;
}
总结
写一半时按一下撤销文章变回半小时小时前的气死我了
这里三角窗是为实现的，用了矩形代替，可自己修改，
可能有bug我没有测试带阻等设计
1491
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
