我們來(lái)實(shí)際操作STM32F4DSP庫(kù)的IIR濾波器。 與IIR濾波器函數(shù)有關(guān)的源文件如下圖所示: STM32F4DSP庫(kù)中采用biquad作為一個(gè)單元。一個(gè)biquad是2階,n個(gè)biquad串聯(lián)之后就是n階濾波器。 基本的單元結(jié)構(gòu)如下所示: 我們可以求出一個(gè)biquad的差分函數(shù)形式是: y[n]=b0*x[n]+b1*x[n-1]+b2*x[n-2]-a1*y[n-1]-a2*y[n-2] Matlab里的計(jì)算就是按照上面的式子計(jì)算的,但是STM32F4DSP庫(kù)里的系數(shù)a1,a2是取反的。 下面就介紹如何使用MATLAB設(shè)計(jì)IIR濾波器。 舉個(gè)例子,我們要設(shè)計(jì)一個(gè)采樣率為1kHz,4階,截止頻率100Hz的巴特沃斯濾波器。 首先和設(shè)計(jì)FIR濾波器一樣,首先在MATLAB命令窗口輸入fdatool,調(diào)出濾波器設(shè)計(jì)窗口 按照方框所示設(shè)置好參數(shù): 點(diǎn)擊DesignFilter: 注意紅框里面是直接II型,我們要把他改為直接I型。 點(diǎn)擊Edit->ConvertStructure,選擇I型: 轉(zhuǎn)化好后,點(diǎn)擊File-Export, 第一項(xiàng)選擇CoefficientFile(ASCII): 第一項(xiàng)選擇好以后,第二項(xiàng)選擇Decimal: 點(diǎn)擊Export,保存后生成如下文件: 系數(shù)對(duì)應(yīng)如下: 1211-1.32091343081942640.63273879288527657 b0b1b2a0a1a2 1211-1.04859957636261170.29614035756166951 b0b1b2a0a1a2 實(shí)際使用ARM官方的IIR函數(shù)調(diào)用的時(shí)候要將a1和a2取反。把a0去掉 ScaleValues表示每個(gè)biquad的增益系數(shù)。所以最后用STM32計(jì)算后,要乘以這兩個(gè)系數(shù)。 設(shè)計(jì)濾波器系數(shù)之后,我們來(lái)看STM32的IIR濾波器函數(shù): 主要介紹下arm_biquad_cascade_df1_f32 函數(shù)定義如下: voidarm_biquad_cascade_df1_f32( constarm_biquad_casd_df1_inst_f32*S, float32_t*pSrc, float32_t*pDst, uint32_tblockSize) 參數(shù): *Spointstoaninstanceofthefloating-pointBiquadcascadestructure. *pSrcpointstotheblockofinputdata. *pDstpointstotheblockofoutputdata. blockSizenumberofsamplestoprocesspercall. 介紹下結(jié)構(gòu)體arm_biquad_casd_df1_inst_f32 typedefstruct { //<numberof2ndorderstagesinthefilter.Overallorderis2*numStages. uint32_tnumStages; //<Pointstothearrayofstatecoefficients.Thearrayisoflength4*numStages.float32_t*pState; //<Pointstothearrayofcoefficients.Thearrayisoflength5*numStages. float32_t*pCoeffs; }arm_biquad_casd_df1_inst_f32; 注意下:pState指向的數(shù)組長(zhǎng)度是4倍numStages長(zhǎng)度 pCoeffs指向的數(shù)組長(zhǎng)度是5倍numStages長(zhǎng)度,a0默認(rèn)為1,不需要放入 numStages表示biquad個(gè)數(shù); 好,接下來(lái)我們就可以使用這個(gè)函數(shù)了 - #define numStages 2
- #define TEST_LENGTH_SAMPLES 1024
- float32_t testInput_f32[TEST_LENGTH_SAMPLES];
- float32_t testOutput[TEST_LENGTH_SAMPLES];
- float32_t IIRStateF32[4*numStages];
- const float32_t IIRCoeffs32LP[5*numStages] =
- {
- 1.0f, 2.0f, 1.0f, 1.3209134308194264f, -0.63273879288527657f,
- 1.0f, 2.0f, 1.0f, 1.0485995763626117f, -0.29614035756166951f
- };
- void arm_iir_f32_lp(void)
- {
- uint32_t i;
- arm_biquad_casd_df1_inst_f32 S;
- float32_t ScaleValue;
- for(i=0;i<TEST_LENGTH_SAMPLES;i++)
- {
- testInput_f32[i]=1.2f*arm_sin_f32(2*PI*50*i/1000)+arm_sin_f32(2*PI*250*i/1000)+1;
-
- printf("%frn", testInput_f32[i]);
- }
- arm_biquad_cascade_df1_init_f32(&S, numStages, (float32_t *)&IIRCoeffs32LP[0], (float32_t *)&IIRStateF32[0]);
- arm_biquad_cascade_df1_f32(&S, testInput_f32, testOutput, TEST_LENGTH_SAMPLES);
- ScaleValue = 0.077956340516462552f * 0.061885195299764481f;
- for(i=0; i<TEST_LENGTH_SAMPLES; i++)
- {
- printf("%frn", testOutput[i]*ScaleValue);
- }
- }
復(fù)制代碼
把原始信號(hào)和過(guò)濾后信號(hào)打印出來(lái),導(dǎo)入到matlab,用下面程序處理: Fs=1000; N=1024; n=0:1:N-1; f=Fs*n/N; t=0:1/Fs:(N-1)/Fs; subplot(2,2,1); plot(t(1:150),data1(1:150)); xlabel('時(shí)間/s'); ylabel('幅度/v'); title('原始信號(hào)波形圖'); h1=fft(data1,N); subplot(2,2,2); plot(t(1:150),data2(1:150)); xlabel('時(shí)間/s'); ylabel('幅度/v'); title('過(guò)濾后信號(hào)波形圖'); subplot(2,2,3); plot(f,abs(h1)); xlabel('頻率/Hz'); ylabel('幅度'); title('原始信號(hào)頻譜圖'); subplot(2,2,4); h2=fft(data2,N); plot(f,abs(h2)); xlabel('頻率/Hz'); ylabel('幅度'); title('過(guò)濾后信號(hào)頻譜圖'); 運(yùn)行結(jié)果:
可以看出STM32的IIR濾波器的計(jì)算結(jié)果還是令人滿意的。
全部資料下載地址:
|