滤波在数学上是如何实现的

在单片机进行数据采集时，会遇到数据的随机误差，随机误差是由随机干扰引起的，其特点是在相同条件下测量同一量时，其大小和符号会现无规则的变化而无法预测，但多次测量的结果符合统计规律。为克服随机干扰引起的误差，硬件上可采用滤波技术，软件上可采用软件算法实现数字滤波。滤波算法往往是系统测控算法的一个重要组成部分，实时性很强。

采用数字滤波算法克服随机干扰的误差具有以下优点：

1、数字滤波无需其他的硬件成本，只用一个计算过程，可靠性高，不存在阻抗匹配问题。尤其是数字滤波可以对频率很低的信号进行滤波，这是模拟滤波器做不到的。

2、数字滤波使用软件算法实现，多输入通道可共用一个滤波程序，降低系统开支。

3、只要适当改变滤波器的滤波程序或运算，就能方便地改变其滤波特性，这对于滤除低频干扰和随机信号会有较大的效果。

4、在单片机系统中常用的滤波算法有限幅滤波法、中值滤波法、算术平均滤波法、加权平均滤波法、滑动平均滤波等。

(1)限幅滤波算法

该运算的过程中将两次相邻的采样相减，求出其增量，然后将增量的绝对值，与两次采样允许的最大差值A进行比较。A的大小由被测对象的具体情况而定，如果小于或等于允许的最大差值，则本次采样有效;否则取上次采样值作为本次数据的样本。

算法的程序代码如下：

#defineA //允许的最大差值

chardata; //上一次的数据

char filter()

{

chardatanew; //新数据变量

datanew=get_data(); //获得新数据变量

if((datanew-data)>A||(data-datanew>A))

return data;

else

returndatanew;

}

说明：限幅滤波法主要用于处理变化较为缓慢的数据，如温度、物体的位置等。使用时，关键要选取合适的门限制A。通常这可由经验数据获得，必要时可通过实验得到。

(2)中值滤波算法

该运算的过程是对某一参数连续采样N次(N一般为奇数)，然后把N次采样的值按从小到大排列，再取中间值作为本次采样值，整个过程实际上是一个序列排序的过程。

算法的程序代码如下：

#define N11 //定义获得的数据个数

char filter()

{

charvalue_buff[N]; //定义存储数据的数组

char count,i,j,temp;

for(count=0;count

{

value_buf[count]=get_data();

delay(); //如果采集数据比较慢，那么就需要延时或中断

}

for(j=0;j

{

for(value_buff[i]>value_buff[i+1]

{

temp=value_buff[i];

value_buff[i]=value_buff[i+1];

value_buff[i+1]=temp;

}

}

returnvalue_buff[(N-1)/2];

}

说明：中值滤波比较适用于去掉由偶然因素引起的波动和采样器不稳定而引起的脉动干扰。若被测量值变化比较慢，采用中值滤波法效果会比较好，但如果数据变化比较快，则不宜采用此方法。

(3)算术平均滤波算法

该算法的基本原理很简单，就是连续取N次采样值后进行算术平均。

算法的程序代码如下：

char filter()

{

int sum=0;

for(count=0;count

{

sum+=get_data();

delay():

}

return (char)(sum/N);

}

说明：算术平均滤波算法适用于对具有随机干扰的信号进行滤波。这种信号的特点是有一个平均值，信号在某一数值附近上下波动。信号的平均平滑程度完全到决于N值。当N较大时，平滑度高，灵敏度低;当N较小时，平滑度低，但灵敏度高。为了方便求平均值，N一般取4、8、16、32之类的2的整数幂，以便在程序中用移位 *** 作来代替除法。

(4)加权平均滤波算法

由于前面所说的“算术平均滤波算法”存在平滑度和灵敏度之间的矛盾。为了协调平滑度和灵敏度之间的关系，可采用加权平均滤波。它的原理是对连续N次采样值分别乘上不同的加权系数之后再求累加，加权系数一般先小后大，以突出后面若干采样的效果，加强系统对参数变化趋势的认识。各个加权系数均小于1的小数，且满足总和等于1的结束条件。这样加权运算之后的累加和即为有效采样值。其中加权平均数字滤波的数学模型是：

式中：D为N个采样值的加权平均值：XN-i为第N-i次采样值;N为采样次数;Ci为加权系数。加权系数Ci体现了各种采样值在平均值中所占的比例。一般来说采样次数越靠后，取的比例越大，这样可增加新采样在平均值中所占的比重。加权平均值滤波法可突出一部分信号抵制另一部分信号，以提高采样值变化的灵敏度。

样例程序代码如下：

char codejq[N]={1,2,3,4,5,6,7,8,9,10,11,12}; //code数组为加权系数表，存在程序存储区

char codesum_jq=1+2+3+4+5+6+7+8+9+10+11+12;

char filter()

{

char count;

char value_buff[N];

int sum=0;

for(count=0;count

{

value_buff[count]=get_data();

delay();

}

for(count=0;count

sum+=value_buff[count]jq[count];

return(char)(sum/sum_jq);

}

(5)滑动平均滤波算法

以上介绍和各种平均滤波算法有一个共同点，即每获取一个有效采样值必须连续进行若干次采样，当采速度慢时，系统的实时得不到保证。这里介绍的滑动平均滤波算法只采样一次，将一次采样值和过去的若干次采样值一起求平均，得到的有效采样值即可投入使用。如果取N个采样值求平均，存储区中必须开辟N个数据的暂存区。每新采集一个数据便存入暂存区中，同时去掉一个最老数据，保存这N个数据始终是最新更新的数据。采用环型队列结构可以方便地实现这种数据存放方式。

程序代码如下：

char value_buff[N];

char i=0;

char filter()

{

char count;

int sum=0;

value_buff[i++]=get_data();

if(i==N)

i=0;

for(count=0;count

sum=value_buff[count];

return (char)(sum/N);

}

(6)低通滤波

将普通硬件RC低通滤波器的微分方程用差分方程来表求，变可以采用软件算法来模拟硬件滤波的功能，经推导，低通滤波算法如下：

Yn=a Xn+(1-a) Yn-1

式中 Xn——本次采样值

Yn-1——上次的滤波输出值;

，a——滤波系数，其值通常远小于1;

Yn——本次滤波的输出值。

由上式可以看出，本次滤波的输出值主要取决于上次滤波的输出值(注意不是上次的采样值，这和加权平均滤波是有本质区别的)，本次采样值对滤波输出的贡献是比较小的，但多少有些修正作用，这种算法便模拟了具体有教大惯性的低通滤波器功能。滤波算法的截止频率可用以下式计算：

fL=a/2Pit pi为圆周率314…

式中 a——滤波系数;

， t——采样间隔时间;

例如：当t=05s(即每秒2次)，a=1/32时;

fL=(1/32)/(231405)=001Hz

当目标参数为变化很慢的物理量时，这是很有效的。另外一方面，它不能滤除高于1/2采样频率的干搅信号，本例中采样频率为2Hz，故对1Hz以上的干搅信号应采用其他方式滤除，

低通滤波算法程序于加权平均滤波相似，但加权系数只有两个：a和1-a。为计算方便，a取一整数，1-a用256-a，来代替，计算结果舍去最低字节即可，因为只有两项，a和1-a，均以立即数的形式编入程序中，不另外设表格。虽然采样值为单元字节(8位A/D)。为保证运算精度，滤波输出值用双字节表示，其中一个字节整数，一字节小数，否则有可能因为每次舍去尾数而使输出不会变化。

设Yn-1存放在30H(整数)和31H(小数)两单元中，Yn存放在32H(整数)和33H(小数)中。滤波程序如下：

虽千万里，吾往矣。

冲击响应不变法函数

[bz,az]=impinvar(b,a,Fs)

[bz,az]=impinvar(b,a)

例如：取采样频率f=1KHz,用双线性变换法设计五阶Butterworth低通数字滤波器，绘出模拟滤波器与数字滤波器的幅频与相频特性，MATLAB程序如下：

[z,p,k]=buttap(5) ;% 设计五阶Butterworth低通模拟滤波器原型

[zd,pd,kd]=bilinear(z,p,k,1000);%双线性变换得到低通数字滤波器

[b,a]=zp2tf(zd,pd,kd);%滤波器类型转换

w=128;

freqs(b,a,w)

figure;

freqz(b,a,w)

实例：1、设带通滤波器的滤波器中心频率为W0=2KHz,带宽为BW=100Hz, 取采样频率f=10kHZ,用脉冲相应不变法设计，设计五阶带通Butterworth数字滤波器,绘出数字滤波器的频谱特性

[z,p,k]=buttap(5);

[b,a]=zp2tf(z,p,k);

w=128;

w0=2000;

[bt,at]=lp2bp(b,a,w0,10000);

[bz,az]=impinvar(b,a,w);

freqz(bt,at,w)

2、直接设计五阶butterworth带通滤波器，绘出频谱图。（高端与低端截止频率分别为02和09）

figure;

w=[02,09];

[b,a]=butter(5,w);

freqz(b,a)

3、设高通截止频率为w0=10000Hz, 取采样频率f=20000,用双线性变换法设计六阶高通Butterworth数字滤波器,绘出数字滤波器的频谱特性

[z,p,k]=besselap(6);

[b,a]=zp2tf(z,p,k);

w=128;

w0=10000;

[bt,at]=lp2hp(b,a,w0);

[bz,az]=bilinear(bt,at,20000);

freqz(bz,az,w)

4、取采样频率f=100Hz,用双线性变换法设计五阶Butterworth低通数字滤波器，绘出模拟滤波器与数字滤波器的幅频与相频特性

[z,p,k]=buttap(5);

[zd,pd,kd]=bilinear(z,p,k,100);

[b,a]=zp2tf(zd,pd,kd);

w=128;

freqs(b,a,w);

figure;

freqz(b,a,w)

具体的得根据情况自己确定

老师给出的程序，大部分，都不用动，

仅仅改动几个地址，即可符合题目要求。

把程序中的：

70H，改为 60H；

71H，改为 61H；

81H，改为 71H；

82H，改为 72H。

即可。

完整的程序

%写上标题

%设计低通滤波器：

[N,Wc]=buttord()

%估算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wc

[a,b]=butter(N,Wc); %设计Butterworth低通滤波器

[h,f]=freqz(); %求数字低通滤波器的频率响应

figure(2); % 打开窗口2

subplot(221); %图形显示分割窗口

plot(f,abs(h)); %绘制Butterworth低通滤波器的幅频响应图

title(巴氏低通滤波器'');

grid; %绘制带网格的图像

sf=filter(a,b,s); %叠加函数S经过低通滤波器以后的新函数

subplot(222);

plot(t,sf); %绘制叠加函数S经过低通滤波器以后的时域图形

xlabel('时间 (seconds)');

ylabel('时间按幅度');

SF=fft(sf,256); %对叠加函数S经过低通滤波器以后的新函数进行256点的基—2快速傅立叶变换

w= %新信号角频率

subplot(223);

plot()); %绘制叠加函数S经过低通滤波器以后的频谱图

title('低通滤波后的频谱图');

%设计高通滤波器

[N,Wc]=buttord()

%估算得到Butterworth高通滤波器的最小阶数N和3dB截止频率Wc

[a,b]=butter(N,Wc,'high'); %设计Butterworth高通滤波器

[h,f]=freqz(); %求数字高通滤波器的频率响应

figure(3);

subplot(221);

plot()); %绘制Butterworth高通滤波器的幅频响应图

title('巴氏高通滤波器');

grid; %绘制带网格的图像

sf=filter(); %叠加函数S经过高通滤波器以后的新函数

subplot(222);

plot(t,sf); ;%绘制叠加函数S经过高通滤波器以后的时域图形

xlabel('Time(seconds)');

ylabel('Time waveform');

w; %新信号角频率

subplot(223);

plot()); %绘制叠加函数S经过高通滤波器以后的频谱图

title('高通滤波后的频谱图');

%设计带通滤波器

[N,Wc]=buttord([)

%估算得到Butterworth带通滤波器的最小阶数N和3dB截止频率Wc

[a,b]=butter(N,Wc); %设计Butterworth带通滤波器

[h,f]=freqz(); %求数字带通滤波器的频率响应

figure(4);

subplot(221);

plot(f,abs(h)); %绘制Butterworth带通滤波器的幅频响应图

title('butter bandpass filter');

grid; %绘制带网格的图像

sf=filter(a,b,s); %叠加函数S经过带通滤波器以后的新函数

subplot(222);

plot(t,sf); %绘制叠加函数S经过带通滤波器以后的时域图形

xlabel('Time(seconds)');

ylabel('Time waveform');

SF=fft(); %对叠加函数S经过带通滤波器以后的新函数进行256点的基—2快速傅立叶变换

w=( %新信号角频率

subplot(223);

plot(')); %绘制叠加函数S经过带通滤波器以后的频谱图

title('带通滤波后的频谱图');

I=imread('lenabmp');

inf=imfinfo('lenabmp')

figure,imshow(I)

X=grayslice(I,64);

imshow(X,pink(64))

load trees

figure,image(10,10,X)

imwrite(X,map,'treesbmp');

imfinfo('treesbmp')

figure,imshow(X,map)

BW=im2bw(X,map,06);

figure,imshow(BW)

I=imread('lenabmp');

inf=imfinfo('lenabmp')

figure,imshow(I)

X=grayslice(I,64);

figure,imshow(X,pink(64))

A=imread('lenabmp');

imshow(A)

B=fftshift(fft2(A));

figure;

imshow(log(abs(B)),[8,10])

clc;

clear all;

I=imread('lenatif');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %用中值滤波，多维滤波，使用中心为-4,-8的拉普

% %拉斯滤波器，高斯低通滤波，拉普拉斯滤波器进行滤波处理

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure;%figure1

subplot(2,2,1);

imshow(I);

title('原始图像');

J=imnoise(I,'salt & pepper',004);%加椒盐噪声

title('加椒盐噪声');

subplot(2,2,2);

imshow(J);

K=medfilt2(J,[4,4])%进行中值滤波;

subplot(2,2,3);

imshow(K);

title('进行中值滤波');

h=ones(3,3)/9;%多维滤波

w=h;

g=imfilter(I,w,'conv','replicate')

subplot(2,2,4);

imshow(g);

title('多维滤波');

%使用中心为-4,-8的拉普拉斯滤波器,

w4=[1 1 1;

1 -4 1;

1 1 1];

w8=[1 1 1;

1 -8 1;

1 1 1];

f=im2double(I);

g4=f-imfilter(f,w4,'replicate');

g8=f-imfilter(f,w8,'replicate');

imshow(f);

figure;%figure2

subplot(1,2,1);

imshow(g4);

title('中心为-4的拉普拉斯滤波');

subplot(1,2,2);

imshow(g8);

title('中心为-8的拉普拉斯滤波');

h3=fspecial('gaussian',[3,3],05);%高斯低通滤波

figure;%figure3

B4=filter2(h3,I);

subplot(1,2,1);

imshow(B4,[ ]);

title('高斯低通滤波');

h4=fspecial('laplacian',0);%使用拉普拉斯滤波器

B5=filter2(h4,I);

subplot(1,2,2);

imshow(B5,[ ]);

title('拉普拉斯滤波器');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %从空域的角度进行亮度变换

% %把灰度等级是10-100的变化到10-255

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure;%figure4

subplot(2,2,1);

imshow(I);

title('原始图像');

J2=imadjust(I,[],[],05);% 增强低灰度级

subplot(2,2,2);

imshow(J2);

title('增强低灰度级');

J3=imadjust(I,[ ],[ ],2);%增强高灰度级

subplot(2,2,3);

imshow(J3);

title('增强高灰度级');

a1=100/255;%把灰度等级是10-100的变化到10-255

a2=255/255;

a3=10/255;

J2=imadjust(I,[a3,a1],[a3,a2],[]);

subplot(2,2,4);

imshow(J2);

title('把灰度等级是10-100的变化到10-255');

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %从频域的角度进行亮度变换

% %fft2

% %由于能量主要集中在低频部分

% %所以对低频进行处理可以得到理想的效果

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I=imread('lenatif');

up=05;%设置处理频率上限

down=009%%设置处理频率下限

figure;%figure5

subplot(421);

imshow(I);

title('原始图像');

f=double(I);

subplot(4,2,3);

imshow(log(abs(f)),[]);

title('unit8转化为double');

g=fft2(f);

g=fftshift(g);

subplot(4,2,5);

imshow(log(abs(g)),[]);

title('FFT2变化后的图像');

[M,N]=size(g);% 转换数据矩阵

y1=max(max(abs(g)));%求出最大频率

y2=min(min(abs(g)));%%求出最小频率

y3=(y1-y2)up+y2;%设置滤波上限

y4=(y1-y2)down+y2;%%设置滤波下限

for i=1:M

for j=1:N

if (abs(g(i,j))<y4)

g(i,j)=g(i,j)^11;%对低频部分进行灰度增强

end

end

end

result=ifftshift(g);

J2=ifft2(result);

J3=uint8(abs(J2));

subplot(427);

imshow(J3,[ ]);

title('频域处理后的图像');

subplot(422)

imhist(I,64);

subplot(424)

imhist(f,64);

subplot(426)

imhist(g,64);

subplot(428)

imhist(J3,64);

clear

N=200;%取200个数

w(1)=0;

w=randn(1,N);%产生一个1×N的行向量，第一个数为0，w为过程噪声（其和后边的v在卡尔曼理论里均为高斯白噪声）

x(1)=0;%状态x初始值

a=1;%a为状态转移阵，此程序简单起见取1

for k=2:N

x(k)=ax(k-1)+w(k-1); %系统状态方程，k时刻的状态等于k-1时刻状态乘以状态转移阵加噪声（此处忽略了系统的控制量）

end

V=randn(1,N);%测量噪声

q1=std(V);

Rvv=q1^2;

q2=std(x);

Rxx=q2^2; %此方程未用到Rxx

q3=std(w);

Rww=q3^2; %Rvv、Rww分别为过程噪声和测量噪声的协方差(此方程只取一组数方差与协方差相同)

c=02;

Y=cx+V;%量测方差，c为量测矩阵，同a简化取为一个数

p(1)=0;%初始最优化估计协方差

s(1)=0;%s(1)表示为初始最优化估计

for t=2:N

p1(t)=a^2p(t-1)+Rww;%p1为一步估计的协方差，此式从t-1时刻最优化估计s的协方差得到t-1时刻到t时刻一步估计的协方差

b(t)=cp1(t)/(c^2p1(t)+Rvv);%b为卡尔曼增益，其意义表示为状态误差的协方差与量测误差的协方差之比(个人见解)

s(t)=as(t-1)+b(t)(Y(t)-acs(t-1));%Y(t)-acs(t-1)称之为新息，是观测值与一步估计得到的观测值之差，此式由上一时刻状态的最优化估计s(t-1)得到当前时刻的最优化估计s(t)

p(t)=p1(t)-cb(t)p1(t);%此式由一步估计的协方差得到此时刻最优化估计的协方差

end

t=1:N;

plot(t,s,'r',t,Y,'g',t,x,'b');%作图，红色为卡尔曼滤波，绿色为量测，蓝色为状态

%整体来说，此卡尔曼程序就是一个循环迭代的过程，给出初始的状态x和协方差p，得到下一时刻的x和p，循环带入可得到一系列的最优的状态估计值，此方法通常用于目标跟踪和定位。

%本人研究方向与此有关，有兴趣可以交流下。

例题：

设计一个低通数字滤波器，给定抽样频率为fs=12000Hz，通带截止频率wp=0Hz，阻带起始频率ws=5000Hz（假设阻带衰减不小于-50dB）。

解答：

由于阻带衰减为50dB，查表，可选海明窗，其阻带最小衰减为53dB，过渡带宽度为66π/N。

MATLAB程序如下：

wp=0;% 频率归一化

ws=50002/12000;

wdel=ws-wp;% 过渡带宽

wn=05(wp+ws);% 近似计算截止频率

N=ceil(66pi/wdel);% 根据过渡带宽度求滤波器阶数

window=hamming(N+1);% 海明窗

b=fir1(N,wn,window);% FIR滤波器设计

freqz(b,1,512);% 查看滤波器幅频及相频特性

参考：

>

你的问题没有说明衰减，我直接在百度上找的答案，使用海明窗的低通滤波器。我运行了下，语句是对的。

以上就是关于滤波在数学上是如何实现的全部的内容，包括:滤波在数学上是如何实现的、求 MATLAB 巴特沃思 低通滤波器程序。、用汇编语言写去极值平均滤波程序等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！