c语言算法a i-j x y,MATLAB里面的filter和filtfilt的C语言源代码
嗯,算法非常简单,就是网上搜不到C代码实现。filter是个很万能的数字滤波器函数,只要有滤波器的差分方程系数,IIR呀FIR呀都能通过它实现。在MATLAB里面,filter最常用的格式是这两个:
[y,zf] = filter(b,a,X)
[y,zf] = filter(b,a,X,zi)
其中b和a就是差分方程的系数,X是输入向量,zi是“初始状态”。可能这么说明还是不很清晰,那么请看图(注意,a(1)为1,这个可以通过差分方程所有系数除以a(1)所得):
嗯,这样子很轻松地就能把各个y值给算出来了,哦注意上面式子里面的n是“向量a或者b中最长的长度”,在实际编程的时候,如果a和b长度不一样,短者显然是要用0补齐的。对于那个初始状态zi,忽略的时候,比如(顺便提醒一句MATLAB的下标从1开始)
y(1)=b(1)x(1)
y(2)=b(1)x(2)+Z1(1)= b(1)x(2) + b(2)x(1) - a(2)y(1)
不忽略的时候
y(1)=b(1)x(1) + Z1(0)
y(2)=b(1)x(2)+Z1(1)= b(1)x(2) + b(2)x(1) - a(2)y(1) + Z2(0)
因为实际程序中自己定义的东西比较多(=,=|||这也是没办法的事情不是),而filtfilt这个超级无敌的“零相移滤波函数”更是复杂到稍微调用了一下自己写的矩阵运算函数,所以代码全部贴上来实在是太乱。等这次工程做完会大概地把代码打包发一下。现在就先贴点代码片段好了……但愿我的代码风格没那么难懂…… #include "filter.h"
#include //Transposed Direct-Form II Realization
//initial conditions: zi, length==nfilt-1. ignore when zi==NULL
#ifndef EPS
#define EPS 0.000001
#endif
void
filter(const double* x, double* y, int xlen, double* a, double* b, int nfilt, double* zi){
double tmp;
int i,j;
//normalization
if( (*a-1.0>EPS) || (*a-1.0 tmp=*a;
for(i=0;ib[i]/=tmp;
a[i]/=tmp;
}
}
memset(y,0,xlen*sizeof(double));
a[0]=0.0;
for(i=0;ifor(j=0;i>=j&&jy[i] += (b[j]*x[i-j]-a[j]*y[i-j]);
}
if(zi&&i}
a[0]=1.0;
}
OK,接下来是神奇的零相移滤波filtfilt,操作上不复杂,原理上小小有点复杂。它主要是先找到一个“合理的”初始状态Zi,使得无论先从反向滤波还是先从正向滤波结果一致,然后filter一次,逆向,再filter一次,再逆向,就是结果了。这里面包括3个要点:
1. Zi的确定。
2. 边缘效应的消除。
3. 正反向滤波的数学原理。
对于要点1,可以参阅Fredrik Gustafsson的论文Determining the Initial States in Forward-backward Filtering的数学证明。要点2,filtfilt对数据两头做了镜像拓延,最后滤完波再截掉头尾。要点3,这个貌似不大难推()
#include #include "filter.h"
#include "mulMat.h"
#include "invMat.h"
int
filtfilt(const double* x, double* y, int xlen, double* a, double* b, int nfilt){
int nfact;
int tlen; //length of tx
int i;
double *tx,*tx1,*p,*t,*end;
double *sp,*tvec,*zi;
double tmp,tmp1;
nfact=nfilt-1; //3*nfact: length of edge transients
if(xlen<=3*nfact || nfilt<2) return -1; //too short input x or a,b
//Extrapolate beginning and end of data sequence using a "reflection
//method". Slopes of original and extrapolated sequences match at
//the end points.
//This reduces end effects.
tlen=6*nfact+xlen;
tx=malloc(tlen*sizeof(double));
tx1=malloc(tlen*sizeof(double));
sp=malloc( sizeof(double) * nfact * nfact );
tvec=malloc( sizeof(double) * nfact );
zi=malloc( sizeof(double) * nfact );
if( !tx || !tx1 || !sp || !tvec || !zi ){
free(tx);
free(tx1);
free(sp);
free(tvec);
free(zi);
return 1;
}
tmp=x[0];
for(p=x+3*nfact,t=tx;p>x;--p,++t) *t=2.0*tmp-*p;
for(end=x+xlen;ptmp=x[xlen-1];
for(end=tx+tlen,p-=2;t//now tx is ok.
end = sp + nfact*nfact;
p=sp;
while(psp[0]=1.0+a[1];
for(i=1;isp[i*nfact]=a[i+1];
sp[i*nfact+i]=1.0L;
sp[(i-1)*nfact+i]=-1.0L;
}
for(i=0;itvec[i]=b[i+1]-a[i+1]*b[0];
}
if(invMat(sp,nfact)){
free(zi);
zi=NULL;
}
else{
mulMat(sp,tvec,zi,nfact,nfact,1);
}//zi is ok
free(sp);free(tvec);
//filtering tx, save it in tx1
tmp1=tx[0];
if(zi)
for( p=zi,end=zi+nfact; pfilter(tx,tx1,tlen,a,b,nfilt,zi);
//reverse tx1
for( p=tx1,end=tx1+tlen-1; ptmp = *p;
*p = *end;
*end = tmp;
}
//filter again
tmp1 = (*tx1)/tmp1;
if(zi)
for( p=zi,end=zi+nfact; pfilter(tx1,tx,tlen,a,b,nfilt,zi);
//reverse to y
end = y+xlen;
p = tx+3*nfact+xlen-1;
while(y*y++ = *p--;
}
free(zi);
free(tx);
free(tx1);
return 0;
}
与MATLAB对比(MATLAB代码):
x=[1 2 3 4 5 6 7 8];
a=[1 2 3];
b=[4 5 6];
t=filtfilt(b,a,x);
for i=1:8, fprintf(1,'%f\n',t(i)), end;
结果为:
-6731884.250000
7501778.750000
-2757230.250000
-662443.250000
1360955.750000
-686678.250000
4135.750000
227147.750000
这个例子用上面给出的C语言版的filtfilt计算结果完全一致:
c语言算法a i-j x y,MATLAB里面的filter和filtfilt的C语言源代码相关推荐
- MATLAB里面的filter和filtfilt的C语言源代码
MATLAB里面的filter和filtfilt的C语言源代码 嗯,算法非常简单,就是网上搜不到C代码实现.filter是个很万能的数字滤波器函数,只要有滤波器的差分方程系数,IIR呀FIR呀都能通过 ...
- C语言实现:输入一串字符把里面的A、a字符替换成C输出
C语言实现:输入一串字符把里面的A.a字符替换成C输出 #include "stdafx.h" #include<stdio.h> #include<string ...
- MATLAB语言算法实验报告,机械工程实验——matlab实验报告.doc
机械工程实验教学中心 - PAGE 20 - 机械工程实验教学中心 实验指导书 实验名称 基于Matlab的信号处理实验 课程名称 自选综合实验 一.实验目的及要求 实验目的 通过基于Matlab的信 ...
- matlab里面的sul,MATLAB语言在电机控制系统仿真研究中的应用
MATLAB语言在电机控制系统仿真研究中的应用 宋凌锋李立毅程树康 [摘要]简要介绍了MATLAB语言,并把MATLAB语言应用于电机控制系统的仿真研究中,同时以一个具体实例较为深入地对其进行了说明. ...
- c语言单片机求最小公倍数,单片机常用的14个C语言算法,要熟记在心哦!
原标题:单片机常用的14个C语言算法,要熟记在心哦! 算法(Algorithm):计算机解题的基本思想方法和步骤. 算法的描述:是对要解决一个问题或要完成一项任务所采取的方法和步骤的描述,包括需要什么 ...
- 100个经典的C语言算法
100个经典的C算法 C语言的学习要从基础开始,这里是100个经典的算法 题目:古典问题:有一对兔子,从出生后第3个月起每个月都生一对兔子,小兔 子长到第三个月后每个月又生一对兔子,假如兔子都不死,问 ...
- c语言二分法_14个经典C语言算法你就不看一眼?(附详细代码)
今天,给大家讲一讲,单片机常用的14个C语言算法(附详细代码)哟! 一.计数.求和.求阶乘等简单算法 此类问题都要使用循环,要注意根据问题确定循环变量的初值.终值或结束条件,更要注意用来表示计数.和. ...
- matlab图像剖线,一种在等值线图上任意截取剖面的Matlab语言算法
第25卷 第3期 2003年8月 物探化探计算技术 V ol 125N o .3 A ug .2003COM PU T I N G T ECHN I Q U ES FOR GEO PH YS I CA ...
- c语言成绩存储的算法思想,C语言算法总结,非常精辟
<C语言算法总结,非常精辟>由会员分享,可在线阅读,更多相关<C语言算法总结,非常精辟(9页珍藏版)>请在人人文库网上搜索. 1.c语言算法摘要2 .牛顿迭代a的平方根[思想摘 ...
最新文章
- 直接通过OptionalAttribute, DefaultParameterValueAttribute定义缺省参数
- 微信小程序修改样式弹框wx.showModal
- 如何从RxJava升级到RxJava2
- Phaser都不懂,还学什么多线程
- 算法复杂度为O(N) 的排序算法
- Django Rest Framework -解析器
- 持续集成与持续部署宝典Part 2:创建持续集成流水线
- CGLib动态代理原理及实现
- Liferay6.2.1 集成 CAS4.0 实现单点登录与应用系统集成
- 0819 - 要想富,追新不守旧
- 解决Nginx + PHP(FastCGI)遇到的502 Bad Gateway错误[原创]
- django WEB聊天室项目
- 实战服务器虚拟化,企业虚拟化实战Vmware篇PDF影印版(张巍著) 56M
- Sprint周期开发总结
- Python与爬虫有什么关系?
- 计算机出现蓝屏cpu很烫,win7系统电脑蓝屏罪魁祸首CPU超频的解决方法
- 中南大学计算机学院复试2021,中南大学2021年硕士研究生拟录取名单汇总
- CTF Crypto---RSA N不互素
- 【高效程序员系列】3、别碰鼠标——让键盘飞起来
- dedecms 织梦配置 手机 wap 站点,并绑定二级域名