第一代心电滤波算法测试
warmonkey2010/07/20电子技术 IP:江西
aVF,带宽1-25Hz
好像FIR低通并没有起作用
IIR高通及50Hz陷波表现良好

测试中
DSCF6003.jpg
液晶截图
aVF_LCD.bmp

AD采样用TIM1比较匹配中断1触发
采样率好像算错了,于是只好乱改一通,用频率计校正,现在是设定为400Hz
原来设置为定时器20分频,预置数2250,理论上是18M/20/2250=400Hz
结果采样率相差很远,大概是160Hz
 
定时器初始化代码,现在的:
  RCC_APB2PeriphClockCmd(RCC_APB2Periph_TIM1,ENABLE);
  
  TIM_DeInit(TIM1);
  
  //APB2=18MHz 
  TIM_TimeBaseStructure.TIM_Prescaler = 10;//原来是20
  TIM_TimeBaseStructure.TIM_CounterMode = TIM_CounterMode_Up;
  TIM_TimeBaseStructure.TIM_Period = 2043; //原来是2205
  //刚好为400sps
  TIM_TimeBaseStructure.TIM_ClockDivision = TIM_CKD_DIV1;
  TIM_TimeBaseStructure.TIM_RepetitionCounter = 0x0;
  TIM_TimeBaseInit(TIM1,&TIM_TimeBaseStructure);
  
  //TIM1 CH1 in Toggle mode
  TIM_OCInitStructure.TIM_OCMode = TIM_OCMode_Toggle;
  TIM_OCInitStructure.TIM_OutputState = TIM_OutputState_Enable;
  TIM_OCInitStructure.TIM_OutputNState = TIM_OutputNState_Disable;
  TIM_OCInitStructure.TIM_Pulse = 1024;
  TIM_OCInitStructure.TIM_OCPolarity = TIM_OCPolarity_Low;
  TIM_OCInitStructure.TIM_OCNPolarity = TIM_OCNPolarity_Low;
  TIM_OCInitStructure.TIM_OCIdleState = TIM_OCIdleState_Reset;
  TIM_OCInitStructure.TIM_OCNIdleState = TIM_OCIdleState_Reset;
  TIM_OC1Init(TIM1,&TIM_OCInitStructure);

不知道为什么会有这个问题


“敏感部位”代码
#include "../platform.h"//变量类型定义
 
#include "iir_coefs.h"
#include "iir_filter.h"
 
static int32 y[IIR_NSEC][3];
static int32 x[IIR_NSEC+1][3];

 //自带缓存,逐个输入数据即可,返回的是滤波后的结果
//伪代码例子:savetosdcard(iir_filter(adc_get()));
int16 iir_filter(int16 in)
{
      uint16 i;
      volatile int64 yout;
      
      x[0][0] = in;
      
      for(i=0;i<IIR_NSEC;i++)
      {
          x[0] = x[0]*IIR_G/IIR_A[0];
          
          y[0] =(x[0]*IIR_B[0]+x[1]*IIR_B[1]+x[2]*IIR_B[2]-y[1]*IIR_A[1]-y[2]*IIR_A[2])/IIR_A[0];
          
          y[2]=y[1];y[1]=y[0];
          x[2]=x[1];x[1]=x[0];
          
          x[i+1][0] = y[0];
      }
      
      yout = x[IIR_NSEC][0]*IIR_G[IIR_NSEC]/IIR_A[IIR_NSEC-1][0];//×îºóÒ»Œ¶µÄÔöÒæ
      
      if(yout>32767) yout=32767;
      if(yout<-32768) yout=-32768;
      return (int16)yout;     
      
}
 //复位滤波器
void iir_reset(void)
{
    uint16 i,j;
    
    for(i=0;i<IIR_NSEC+1;i++)
    {
       for(j=0;j<3;j++)
       {
          x[j]=0;
       }
    }
    
    for(i=0;i<IIR_NSEC;i++)
    {
       for(j=0;j<3;j++)
       {
          y[j]=0;
       }
    }
}

系数使用说明
FIR就不提了,说下IIR

用fdatool设计好以后,要转换(convert structure)为Direct-form I而且要是由2阶部分(section)串联的(默认就是这样,不用特别设置)
然后把系数四舍五入(quantize),使系数变为整数,一般选择fixed point ,word length =16,
并且确定滤波器稳定,再导出(export ->export to workspace)
还要导出为C文件一次
%SOS Martix与系数的对应关系
%SOS(j,i)
%在第j个Section中:   
%SOS(j,1:3) b0~b2   别名: NUM[j][0]~NUM[j][3]       B[][]
%SOS(j,4)   a0                  别名: DEN[j][0]                       A[][0]
%SOS(j,5:6) a1~a2    别名:DEN[j][1]~DEN[j][2]       A[][]
%G(i)是第i个Section的输入增益
%G的最后一个元素是输出增益
%Direct-Form I IIR Filter 的差分方程
% a0*y = b0*x + b1*x_1 + b2*x_2 - a1*y_1 -a2*y_2;
把C文件中的NUM,DEN数组的奇数项去掉,并改名为B,A数组,保存为iir_coefs.h给iir_filter.c使用
G矩阵要全部乘以A【任意项(每一项都是相同的)】【0】,以转化为整数,然后写入数组

宏IIR_NSEC是Section(部分)数,必须定义

例子:
//iir_coefs.h
#define IIR_NSEC 5
const int16 IIR_B[IIR_NSEC][3] = {
  {
    16384, -32768,  16384 
  },
  {
    16384, -32767,  16384 
  },
  {
    16384, -32766,  16384 
  },
  {
    16384, -32766,  16384 
  },
  {
    16384, -23170,  16384    
  } 
};
 
const int16 IIR_A[IIR_NSEC][3] = {
 
  {
    16384, -31753,  15404 
  },
  {
    16384, -32596,  16221 
  },
  {
    16384, -32717,  16338 
  },
  {
    16384, -32753,  16373 
  },
  {
    16384, -21802,  14449
  }
};
 
const int16 IIR_G[IIR_NSEC+1] = {
    14158,16301,16360,16378,15417,16384
      };
 
来自:电子信息 / 电子技术
12
已屏蔽 原因:{{ notice.reason }}已屏蔽
{{notice.noticeContent}}
~~空空如也
warmonkey 作者
14年7个月前 IP:未同步
238460
占位,明天附上fdatool生成的系数转换为这套代码能使用的头文件的方法。
另外在此感谢各位网友的帮助,尤其是小俊
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
kknd
14年7个月前 IP:未同步
238466
又有进展了。。。
恭喜恭喜。。。 [s:252]
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
托尼史塔克
14年7个月前 IP:未同步
238507
这是,这是,这是干嘛用的?
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小特斯拉圈圈
14年7个月前 IP:未同步
238534
引用第3楼托尼史塔克于2010-07-21 10:00发表的  :
这是,这是,这是干嘛用的?

LZ DIY心电仪?
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小特斯拉圈圈
14年7个月前 IP:未同步
238535
那个液晶分辨率太低了..
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
warmonkey作者
14年7个月前 IP:未同步
238608
320*240的那块被我的贱手弄烂了,
另外滤波器似乎又有振荡的倾向,折腾中
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
小哈_abc
14年7个月前 IP:未同步
238702
诺基亚手机上的拆机或者诺基亚库存余下的液晶屏吧~~~~~~~~

很便宜的说
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
blackiller
14年7个月前 IP:未同步
238707
确实很牛呀
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
三硝基甲苯
14年7个月前 IP:未同步
238712
额。。。这就是LZ大名鼎鼎的心电图机???佩服。。。机芯是不错,液晶屏太差了。。。
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
warmonkey作者
14年7个月前 IP:未同步
238849
今天出了点岔子,底噪不知怎么变大了很多,
估计是程序问题,也可能是模拟电路的故障,示波器不在手边,难办。
IIR高通改为全浮点计算,能直接支持matlab生成的系数,哇哈哈。

附上代码(iir_filter.c):

#include "datatype.h"
#include "iir_filter.h"
#include "iir_coefs.h"

static float y[IIR_NSEC][3];

static float x[IIR_NSEC+1][3];



int16 iir_filter(int16 in)

{

      uint16 i;

      

      x[0][0] = in;

      

      for(i=0;i<IIR_NSEC;i++)

      {

        

          y[i][0] =x[i][0]*IIR_B[i][0]+x[i][1]*IIR_B[i][1]+x[i][2]*IIR_B[i][2]-y[i][1]*IIR_A[i][1]-y[i][2]*IIR_A[i][2];

          y[i][0] /= IIR_A[i][0];

          

          y[i][2]=y[i][1];y[i][1]=y[i][0];

          x[i][2]=x[i][1];x[i][1]=x[i][0];

          

          x[i+1][0] = y[i][0];

      }

      

      if( x[IIR_NSEC][0]>32767)  x[IIR_NSEC][0]=32767;

      if( x[IIR_NSEC][0]<-32768) x[IIR_NSEC][0]=-32768;

      return  ((int16)x[IIR_NSEC][0]);    

      

}


//复位滤波器
void iir_reset(void)

{

    uint16 i,j;

    

    for(i=0;i<IIR_NSEC+1;i++)

    {

       for(j=0;j<3;j++)

       {

          x[i][j]=0;

       }

    }

    

    for(i=0;i<IIR_NSEC;i++)

    {

       for(j=0;j<3;j++)

       {

          y[i][j]=0;

       }

    }

}

系数生成方法
用fdatool设计一个IIR滤波器,然后点击菜单中的Edit->Convert Structure 选择Direct Form I ,SOS,(必须是Direct Form I,II不行)
注意这里千万不要点击Convert to single section 这一项。
然后选择quantize filter,精度选择single precision floating point (单精度浮点)
再点击菜单中的Targets -> generate c header ,选择export as:single precision floating point (单精度浮点)
填写变量名称时,把NUM改成IIR_B,DEN改成IIR_A,其他不用动,保存为iir_coefs.h
最后打开iir_coefs.h把MWSPT_NSEC替换成IIR_NSEC,   NL、DL删除掉,real32_T改成float ,
其中有一个#include项要删掉,那个twwtypes.h没必要使用的.
改完的文件:
#define IIR_NSEC 9
   //原来叫做MWSPT_NSEC


const float IIR_B[IIR_NSEC][3] = {
   //为什么改为float很明显了吧
  {

     0.8641357422,              0,              0

  },

  {

                1,             -2,              1

  },

  {

     0.9949035645,              0,              0

  },

  {

                1,   -1.999938965,              1

  },

  {

     0.9985351563,              0,              0

  },

  {

                1,    -1.99987793,              1

  },

  {

     0.9996337891,              0,              0

  },

  {

                1,    -1.99987793,              1

  },

  {

                1,              0,              0

  }

};



const float IIR_A[IIR_NSEC][3] = {

  {

                1,              0,              0

  },

  {

                1,   -1.938049316,   0.9401855469

  },

  {

                1,              0,              0

  },

  {

                1,   -1.989501953,   0.9900512695

  },

  {

                1,              0,              0

  },

  {

                1,   -1.996887207,   0.9971923828

  },

  {

                1,              0,              0

  },

  {

                1,   -1.999084473,   0.9993286133

  },

  {

                1,              0,              0

  }

};


最后调用iir_filter函数,就可以得到滤波后的结果
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论
冲压发动机
14年7个月前 IP:未同步
239872
LZ真NB!顶!
引用
评论
加载评论中,请稍候...
200字以内,仅用于支线交流,主线讨论请采用回复功能。
折叠评论

想参与大家的讨论?现在就 登录 或者 注册

所属专业
所属分类
上级专业
同级专业
warmonkey
学者 机友
文章
363
回复
8005
学术分
12
2008/10/11注册,3时10分前活动

Cubesat

主体类型:个人
所属领域:无
认证方式:手机号
IP归属地:未同步
插入公式
评论控制
加载中...
文号:{{pid}}
投诉或举报
加载中...
{{tip}}
请选择违规类型:
{{reason.type}}

空空如也

加载中...
详情
详情
推送到专栏从专栏移除
设为匿名取消匿名
查看作者
回复
只看作者
加入收藏取消收藏
收藏
取消收藏
折叠回复
置顶取消置顶
评学术分
鼓励
设为精选取消精选
管理提醒
编辑
通过审核
评论控制
退修或删除
历史版本
违规记录
投诉或举报
加入黑名单移除黑名单
查看IP
{{format('YYYY/MM/DD HH:mm:ss', toc)}}