FIR數字濾波器 C語言頻域實作 matlab設計FIR參數

• 一、FIR實作步驟
• • 1、FIR參數做FFT
  • 2、輸入資料做FFT
  • 3、FIR做複數乘積
  • 4、将複數乘積結果做IFFT輸出
• 二、FIR代碼實作
• • 1、FIR參數的生成
  • 2、FIR系數做FFT計算
  • 3、輸入資料做overlap
  • 4、輸入資料做FFT
  • 5、資料做複數乘積
  • 6、資料做IFFT并輸出
• 三、代碼說明以及待完善
• • 1、說明
  • 2、待完善

本文主要介紹用C語言實作頻域FIR。

一、FIR實作步驟

1、FIR參數做FFT

2、輸入資料做FFT

3、FIR做複數乘積

4、将複數乘積結果做IFFT輸出

二、FIR代碼實作

1、FIR參數的生成

FIR的系數用matlab的FDA工具生成，設計參數：

FS：48000Hz

Fpass:2000Hz

Fstop:2500Hz

濾波器：低通 FIR

FIR階數：511

Density Factor:20

參數設定如下圖：

2、FIR系數做FFT計算

void FIRCoffesCal(void *processor)
{
	FFT_ELEM *fft = (FFT_ELEM *)processor;
	int i = 0,j = 0;
	memset(fft->overlap   ,0,sizeof(float)*(FFT_LENGTH<<1));
	memset(fft->FIR_CoffesIn,0,sizeof(float)*(FFT_LENGTH<<1));
	for(i=0; i< FIR_ORDER; i++)
	{
		fft->FIR_CoffesIn[i<<1]	= FIR_COFFES[i] ;//* fft->win[i]; //實部
	}
	DSPF_sp_fftSPxSP_cn(FFT_LENGTH, fft->FIR_CoffesIn,  fft->tw , fft->FIR_CoffesOut,fft->base, 0,FFT_LENGTH);
 
}
           

3、輸入資料做overlap

for(i=0; i< FFT_LENGTH-FRAME_LEN; i++)
	{
		fft->frameIn[i]   = fft->frameIn[i+FRAME_LEN] ;
	}
	for(i=0; i< FRAME_LEN; i++)
	{
		fft->frameIn[FFT_LENGTH-FRAME_LEN+i]   = x[i] ;
	}
           

4、輸入資料做FFT

for(i=0; i< FFT_LENGTH; i++)
	{
		fft->FFT_InData[i<<1]   = fft->frameIn[i] ;//* fft->win[i]; //實部
		fft->FFT_InData[(i<<1)+1] = 0 ;			   //虛部賦0
	}
	DSPF_sp_fftSPxSP_cn(FFT_LENGTH, fft->FFT_InData,  fft->tw , fft->FFT_outData,fft->base, 0,FFT_LENGTH);
           

5、資料做複數乘積

/*
struct Complex MUL(struct Complex a,struct Complex b)//定義複乘
{
   struct Complex c;
   c.real=a.real*b.real-a.imag*b.imag;
   c.imag=a.real*b.imag+a.imag*b.real;
   return(c);
}
*/
    fft->FIR_InData[0] = fft->FFT_outData[0]*fft->FIR_CoffesOut[0] - fft->FFT_outData[1]*fft->FIR_CoffesOut[1];//8.569;//
	fft->FIR_InData[1] = fft->FFT_outData[0]*fft->FIR_CoffesOut[1] + fft->FFT_outData[1]*fft->FIR_CoffesOut[0];
	fft->FIR_InData[FFT_LENGTH] = fft->FFT_outData[FFT_LENGTH]*fft->FIR_CoffesOut[FFT_LENGTH] - fft->FFT_outData[FFT_LENGTH+1]*fft->FIR_CoffesOut[FFT_LENGTH+1];
	fft->FIR_InData[FFT_LENGTH+1] = fft->FFT_outData[FFT_LENGTH]*fft->FIR_CoffesOut[FFT_LENGTH+1] + fft->FFT_outData[FFT_LENGTH+1]*fft->FIR_CoffesOut[FFT_LENGTH];

	for(i=2;i<FFT_LENGTH;i+=2)
	{
		fft->FIR_InData[i]   = fft->FFT_outData[i]*fft->FIR_CoffesOut[i]   - fft->FFT_outData[i+1]*fft->FIR_CoffesOut[i+1];
		fft->FIR_InData[i+1] = fft->FFT_outData[i]*fft->FIR_CoffesOut[i+1] + fft->FFT_outData[i+1]*fft->FIR_CoffesOut[i];
		fft->FIR_InData[FFT_LENGTH*2-i]   =   fft->FIR_InData[i];
		fft->FIR_InData[FFT_LENGTH*2-i+1] =  -fft->FIR_InData[i+1];
	}
           

6、資料做IFFT并輸出

DSPF_sp_ifftSPxSP_cn(FFT_LENGTH, fft->FIR_InData,  fft->tw , fft->FIR_OutData,fft->base, 0,FFT_LENGTH);

	for(j=0;j<FRAME_LEN;j++)
	{
		y[j] = fft->FIR_OutData[(FFT_LENGTH>>1)+j<<1];//取後半部分
	}
           

三、代碼說明以及待完善

1、說明

代碼實作了頻域FIR，針對overlap有兩種方案，一種是FFT前做，一種是IFFT後做。代碼未貼出，讀寫檔案操作，以
及FFT和IFFT的實作，這個實作網上應該有大把。其中最關鍵的地方就是FIR複數乘積計算。尤其注意資料對稱問題。
           

2、待完善

代碼暫時支援的濾波器階數必須小于FFT的長度，是以FFT長度做了一個特殊處理。後續完成濾波器長度大于FFT長度的實作。