|
|
matlab的fdatool是好東西，不過很多人不知道該怎么使用它生成的C頭文件。
趁著放假有時間，摸索了幾天，終于搞定。
該程序已經用于心電采集實驗
導聯aVF，帶寬1-25Hz
![]()
實驗過程中圖片
file:///C:/Users/WARMON%7E1/AppData/Local/Temp/moz-screenshot-2.png![]()
屏顯
不多說，切入正題
這里有個fdatool設計的IIR高通濾波器，采樣率400Hz時截止頻率1Hz。
設計定型之后，要做些調整。
以下說明中的英文名詞有些可能對不上fdatool界面上的原文，請大家意會吧
第一步：
點擊菜單中的Edit->Convert Structure 選擇Direct Form I ，SOS，（必須是Direct Form I， II不行）
一般情況下，按照默認設置，fdatool設計都是由二階部分串聯組成的。
這種結構的濾波器穩定性比一個section的要好很多，其他方面的性能也好些。
如果不是的話，點擊Convert to second order sections。
這時，濾波器的結構（structure）應該顯示為 Direct Form I，second order sections
第二步：
選擇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
保存好的文件如下：
//一大堆注釋
//然后：
/* General type conversion for MATLAB generated C-code */
#include "tmwtypes.h"
/*
* Expected path to tmwtypes.h
* C:\Program Files\MATLAB\R2010a\extern\include\tmwtypes.h
*/
/*
* Warning - Filter coefficients were truncated to fit specified data type.
* The resulting response may not match generated theoretical response.
* Use the Filter Design & Analysis Tool to design accurate
* single-precision filter coefficients.
*/
#define MWSPT_NSEC 9
const int NL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1 };
const real32_T IIR_B[MWSPT_NSEC][3] = {
{
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 int DL[MWSPT_NSEC] = { 1,3,1,3,1,3,1,3,1 };
const real32_T IIR_A[MWSPT_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_coefs.h把MWSPT_NSEC替換成IIR_NSEC,
NL、DL數組刪除掉，real32_T改成float ，
其中有一個#include "twmtypes.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
}
};
保存文件，然后使用以下代碼進行濾波
這段代碼是根據Direct Form I 2階IIR濾波的差分方程編寫的
a0*y[n] = b0*x[n] + b1*x[n-1] + b2*x[n-2] - a1*y[n-1] -a2*y[n-2];
//iir_filter.c
#include "../platform.h"
//#include "iir_coefs.float.flat.h"
//#include "iir_coefs.float.sharp.h"
#include "iir_coefs_pass@2Hz_stop@0.8Hz.h"
#include "iir_filter.h"
static float y[IIR_NSEC][3];
static float x[IIR_NSEC+1][3];
//IIR_NSEC階直接型II IIR濾波器
//IIR_NSEC個二階biquad串聯
int16 iir_filter(int16 in)
{
uint16 i;
x[0][0] = in;
for(i=0;i
{
// 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];
y[0] = 0;
if(IIR_B[0] == 1) y[0]+=x[0];
else if(IIR_B[0] == -1) y[0]-=x[0];
else if(IIR_B[0] == -2) y[0]=y[0]-x[0]-x[0];
else if(IIR_B[0] == 0);
else y[0] += x[0]*IIR_B[0];
if(IIR_B[1] == 1) y[0]+=x[1];
else if(IIR_B[1] == -1) y[0]-=x[1];
else if(IIR_B[1] == -2) y[0]=y[0]-x[1]-x[1];
else if(IIR_B[1] == 0);
else y[0] += x[1]*IIR_B[1];
if(IIR_B[2] == 1) y[0]+=x[2];
else if(IIR_B[2] == -1) y[0]-=x[2];
else if(IIR_B[2] == -2) y[0]=y[0]-x[2]-x[2];
else if(IIR_B[2] == 0);
else y[0] += x[2]*IIR_B[2];
if(IIR_A[1] == 1) y[0]-=y[1];
else if(IIR_A[1] == -1) y[0]+=y[1];
else if(IIR_A[1] == -2) y[0]=y[0]+y[1]+y[1];
else if(IIR_A[1] == 0);
else y[0] -= y[1]*IIR_A[1];
if(IIR_A[2] == 1) y[0]-=y[2];
else if(IIR_A[2] == -1) y[0]+=y[2];
else if(IIR_A[2] == -2) y[0]=y[0]+y[2]+y[2];
else if(IIR_A[2] == 0);
else y[0] -= y[2]*IIR_A[2];
if(IIR_A[0] != 1) y[0] /= 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];
}
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
{
for(j=0;j<3;j++)
{
x[j]=0;
}
}
for(i=0;i
{
for(j=0;j<3;j++)
{
y[j]=0;
}
}
}
//iir_filter.h
#ifndef _IIR_FILTER_H__
#define _IIR_FILTER_H__
int16 iir_filter(int16 x);
void iir_reset(void);
#endif
使用方法:
首先寫好iir_coefs.h，然后調用iir_filter.c對數據流進行濾波
一個偽代碼例子：
while(運行中)
{
保存到SD卡(iir_filter(讀取ADC采樣值()));
}
這個函數比STM32 DSP庫中的函數要好很多，DSP庫中的2個IIR濾波函數都不能連續處理數據流。
記得在開始濾波之前重置濾波器
iir_reset();
另：還可以用simulink導出代碼的方式解決，不過這個蛋疼的代碼既然弄出來了，就用著吧。。。
|
|
發表于 2011-1-30 19:51:33
收藏
頂
踩