一阶数字高通滤波器
本文的主要内容包含一阶高通滤波器公式的推导和数字算法的实现以及编程和仿真
1 计算公式推导
1.1.2 算法实现及仿真
利用python实现的代码如下：
import numpy as np
# from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
#2pifWc = 2*np.pi*30Tsw = 0.00314 #采样时间
halfdigiW=np.tan(Wc/2*Tsw )b10=1/(halfdigiW+1)
b11=-b10
a10=(1-halfdigiW)/(halfdigiW+1)x=np.linspace(-np.pi,np.pi,2000) #在[-pi,pi]区间上分割正2000个点 可以理解为信号采样时间为 2*pi/2000s
data=np.zeros_like(x) #输入信号 保存被干扰的信号
data1=np.zeros_like(x) #输入信号,保存未被干扰的信号，方便与滤波后的波形进行比较
y1=np.zeros_like(x) #一阶滤波输出
y2=np.zeros_like(x) #二阶滤波输出 陷波：对希望滤除的特定频率有很好的滤除作用
for i in range(len(x)):data[i] =np.sin( 2 * np.pi * x[i])+0.5*np.sin(30* 2 * np.pi * x[i]) #幅值为1频率为 1Hz的低频信号 + 幅值为0.5 频率为30hz的高频信号data1[i]=0.5*np.sin(130* 2 * np.pi * x[i])# y2[i] = b0*data[i]+b1*data[i-1]+a0*y2[i-1]y1[i] = b10 * data[i] + b11 * data[i - 1] + a10 * y1[i - 1]
y2[0] = 0
y2[1] = 0
y1[0] = 0
y1[1] = 0
#绘原始信号 + 滤波后的信号
plt.subplot(2, 1, 1)
plt.plot(x,data ,label='sig+noise')
plt.plot(x,y1, 'r',label='first order HP')#绘制理想信号 + 滤波后的信号
plt.subplot(2, 1, 2)
plt.plot(x,data1 ,label='sig')
plt.plot(x,y1, 'r',label='first order HP')
plt.grid()
plt.legend()
plt.show()以下是在python中仿真的波形图：
图1-2 fc=0.5Hz
图1-3 fc=2Hz
图1-4 fc=5Hz
图1-5 fc=10Hz
图1-6 fc=20Hz
图1-7 fc=30Hz
图1-8 fc=40Hz
图1-9 fc=60Hz
图1-10 fc=159Hz
输入的信号是幅值为1 频率为 1Hz的低频信号 加上 幅值为0.5 频率为30Hz的高频信号，采样时间为0.003s，从图1-2到图1-10的仿真波形可以看出，当fc为0.5Hz时，滤后的波形有微小的衰减作用，但几乎和原波形一致，当逐渐增大截止频率fc，对低频的滤除结果越来越强，高频越来越接近高频本身的波形，当截止频率高于高频频率时，高频本身也会被滤掉。当截止频率大于等于1/2采样频率时，输出的是一条直线。
利用C语言实现的代码如下：
#ifndef _MHPF1W_F_H_
#define _MHPF1W_F_H_#include <stdint.h>
struct MHpf1W_F
{/*初始化*/struct{void (*Set)(struct MHpf1W_F *self, float cutFreq, int samFreq); //设置截止和采样频率void (*VaryCutFreq)(struct MHpf1W_F *self, float cutFreq); //改变截止频率float cutFreq; //截止频率float samFreq; //采样频率} Init;/*采样计算*/struct{int (*In)(struct MHpf1W_F *self, int Xn);int out_y; //输出值} Prd;/*变量 中间变量 系数等，由初始参数 初始化计算得出*/struct{float Ts; //采样周期int a0, b0, b1; //差分系数int Xn_1, Yn_1;} pri;
};
void MHpf1W_F_Create(struct MHpf1W_F *self);
#endif
//创建方式
// struct MHpf1W_F mlp;
// MHpf1W_F_Create(&mlp);
// mlp.Init.Set(&mlp,2, 1000); fc=20Hz fs=1000Hz#include "MHpf1W_F.h"
#include <string.h>
#include "math.h"static const float PI = 3.1415926535897932384626f;
#define MID(a,min,max) (a= (a<min)?min:(a<max)?a:max)
#define Q15M(a,b) ((a*b)>>15)
/******************************************************************************** 函 数 名 : _Update* 函数功能 : 各系数计算，参数更新* 输入参数 : 滤波器对象 struct MHpf1W_F *self* 返 回 值 : 无*******************************************************************************/
static void _Update(struct MHpf1W_F *self)
{float halfdigiW,tgAnaWT ;halfdigiW = PI * self->Init.cutFreq * self->pri.Ts;tgAnaWT = tan(halfdigiW); //ignore the 1/Tself->pri.b0 = 1/(tgAnaWT+1)*32768; //转成Q15格式self->pri.b1 = -self->pri.b0; //转成Q15格式self->pri.a0 =((1-tgAnaWT)/(tgAnaWT+1))*32768; //转成Q15格式self->pri.Xn_1 = 0;self->pri.Yn_1 = 0;
}
/******************************************************************************** 函 数 名 : InitSet* 函数功能 : 初始化* 输入参数 : cutFreq----截至频率* samFreq----计算机采样频率* 返 回 值 : 无*******************************************************************************/
static void InitSet(struct MHpf1W_F *self, float cutFreq, int samFreq)
{self->Init.cutFreq = MID(cutFreq , 0.0f , samFreq*0.5f);; //截止频率self->Init.samFreq = samFreq; //采样频率self->pri.Ts = 1.0f / self->Init.samFreq; //采样周期 1/Ts_Update(self);
}
/******************************************************************************** 函 数 名 : InitVaryCutF* 函数功能 : 改变截止频率* 输入参数 : cutFreq----截至频率* 返 回 值 : 无*******************************************************************************/
static void InitVaryCutF(struct MHpf1W_F *self, float cutFreq)
{self->Init.cutFreq = cutFreq;_Update(self);
}
/******************************************************************************** 函 数 名 : PrdIn* 函数功能 : 本次输出结果计算* 输入参数 : Xn----本次输入值* 返 回 值 : 本次滤波后的值* 计算公式 ：Y(n)=b0*X(n)+b1*X(n-1)+a0*Y(n-1)*******************************************************************************/
static int PrdIn(struct MHpf1W_F *self, int Xn)
{/*Y(n)=b0*X(n)+b1*X(n-1)+b2*X(n-2)+a0*Y(n-1)+a2*Y(n-2)*/self->Prd.out_y = Q15M(self->pri.b0 , Xn ) + \Q15M(self->pri.b1 , self->pri.Xn_1 ) + \Q15M(self->pri.a0 , self->pri.Yn_1 ) ;self->pri.Yn_1 = self->Prd.out_y;self->pri.Xn_1 = Xn;return self->Prd.out_y;
}
/******************************************************************************** 函 数 名 : MHpf1W_F_Create* 函数功能 : 创建对象 初始化* 输入参数 : self对象* 返 回 值 : 无*******************************************************************************/
void MHpf1W_F_Create(struct MHpf1W_F *self)
{memset(self, 0, sizeof(struct MHpf1W_F));self->Init.Set = InitSet;self->Init.VaryCutFreq = InitVaryCutF;self->Prd.In = PrdIn;
}单片机+匿名科创地面站的软件输出波形如下：
图1-11 fc=0.5Hz
图4-12 fc=2Hz
图1-13 fc=5Hz
图1-14 fc=10Hz
图1-15 fc=20Hz
图1-16 fc=30Hz
图1-17 fc=40Hz
图1-18 fc=60Hz
图1-18 fc=100Hz
图1-18 fc=200Hz
图1-18 fc=500Hz
单片机模拟输入的信号是幅值为1000 频率为 1Hz的低频信号加上幅值为500 频率为30Hz的高频信号，采样时间为0.001s，从图1-11到图1-18的波形可以看出，输出变化的规律现象和python仿真的规律和现象是一致的，同样当截止频率大于等于1/2采样频率时，输出的是一条直线。