天天看點

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

1.機關沖擊響應與頻響

        就如同之前所說的一樣,使用下圖所示的機關沖擊響應,所設計的濾波器,是無法實作的。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

         現在,讓我們看看其這個濾波器的頻響。所謂頻響,就是計算其機關沖擊響應的離散時間傅裡葉變換,

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

       我們可以看出,這個濾波器的頻響的計算結果是實數,并沒有虛數部分。也就是,其相位譜一直是0,也就意味着,這個濾波器輸入與輸出之間沒有延遲,這種相位特性稱為0延遲相位特性。

        但是,這個濾波器無法是無法實作的。我們實際計算一下該濾波器的輸入輸出,就可以明白。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

        這個濾波器在計算的過程中,需要過去的值和未來的值。未來的值是不可預測的,是以,這個濾波器無法實作。為了使得這個濾波器可以實作,我們隻好移動其機關沖擊響應,使得其不再需要未來的值。比如,就像下面這樣移動。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

        這樣的話,這個濾波器就可以實作了,我們隻需要記錄其40個過去的輸入值和現在的輸入值。但是,由于移動了其機關沖擊響應,會不會對頻響産生什麼影響呢,我們來看看。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

       為了更好的說明問題,L去代替之前例子中的20。

       移動之後頻響,我們根據上面式子可以得出兩個結論:1,移動不會對幅度譜産生影響。2,,移動會對相位産生一個延遲,這個延遲主要取決于移動的長度,移動的長度越長,延遲越大。但是,這個移動是線性的。

       是以,我們把這個移動的相位特性稱為,線性相位特性。到這裡,我們移動後的,因果的,可實作的濾波器的機關沖擊響應,如下所示。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

2.窗函數實作的FIR濾波器代碼(C語言)

#include <stdio.h>
#include <math.h>
#include <malloc.h>
#include <string.h>


#define   pi         (3.1415926)

/*-------------Win Type----------------*/
#define   Hamming    (1)



double Input_Data[] = 
{
0.000000 , 0.896802 , 1.538842 , 1.760074 ,  1.538842 ,  1.000000 ,  0.363271 , -0.142040 , -0.363271 , -0.278768,
0.000000 , 0.278768 , 0.363271 , 0.142020 , -0.363271 , -1.000000 , -1.538842 , -1.760074 , -1.538842 , -0.896802,
0.000000 , 0.896802 , 1.538842 , 1.760074 ,  1.538842 ,  1.000000 ,  0.363271 , -0.142040 , -0.363271 , -0.278768,
0.000000 , 0.278768 , 0.363271 , 0.142020 , -0.363271 , -1.000000 , -1.538842 , -1.760074 , -1.538842 , -0.896802,
0.000000 , 0.896802 , 1.538842 , 1.760074 ,  1.538842 ,  1.000000 ,  0.363271 , -0.142040 , -0.363271 , -0.278768,
0.000000 , 0.278768 , 0.363271 , 0.142020 , -0.363271 , -1.000000 , -1.538842 , -1.760074 , -1.538842 , -0.896802,
0.000000 , 0.896802 , 1.538842 , 1.760074 ,  1.538842 ,  1.000000 ,  0.363271 , -0.142040 , -0.363271 , -0.278768,
0.000000 , 0.278768 , 0.363271 , 0.142020 , -0.363271 , -1.000000 , -1.538842 , -1.760074 , -1.538842 , -0.896802,
0.000000 , 55
};




double sinc(double n)
{
    if(0==n) return (double)1;
    else return (double)(sin(pi*n)/(pi*n));
}
 
int Unit_Impulse_Response(int N,double w_c,
                          int Win_Type,
                          double *Output_Data)
{
    signed int Count = 0; 
  
    for(Count = -(N-1)/2;Count <= (N-1)/2;Count++)
    {
    	*(Output_Data+Count+((N-1)/2)) = (w_c/pi)*sinc((w_c/pi)*(double)(Count));
    	//printf("%d %lf  ",Count+((N-1)/2)+1,*(Output_Data+Count+((N-1)/2)));
        //if(Count%4 == 0) printf("\n");
    }	
	
	
    switch (Win_Type)
    {
    	
    	case Hamming:   printf("Hamming \n");
    	                for(Count = -(N-1)/2;Count <= (N-1)/2;Count++)
    			{
    			    *(Output_Data+Count+((N-1)/2)) *= (0.54 + 0.46 * cos((2*pi*Count)/(N-1)));
    			    //printf("%d %lf  ",Count+((N-1)/2)+1,*(Output_Data+Count+((N-1)/2)));
        		    //if(((Count+1)%5 == 0)&&(Count != -20)) printf("\n");
   			} 
    	                break;
    	                
    	                
    	default:        printf("default Hamming \n");
    	                for(Count = -(N-1)/2;Count <= (N-1)/2;Count++)
    			{
    			    *(Output_Data+Count+((N-1)/2)) *= (0.54 + 0.46 * cos((2*pi*Count)/(N-1)));
    			    //printf("%d %lf  ",Count+((N-1)/2)+1,*(Output_Data+Count+((N-1)/2)));
        		    //if(((Count+1)%5 == 0)&&(Count != -20)) printf("\n");
   			} 
    	
                        break;
    }
    
    return (int)1;
}


void Save_Input_Date (double Scand,
                      int    Depth,
                      double *Input_Data)
{
    int Count;
  
    for(Count = 0 ; Count < Depth-1 ; Count++)
    {
    	*(Input_Data + Count) = *(Input_Data + Count + 1);
    }
    
    *(Input_Data + Depth-1) = Scand;
}



double Real_Time_FIR_Filter(double *b,
                            int     b_Lenth,
                            double *Input_Data)
{    
    int Count;
    double Output_Data = 0;
    
    Input_Data += b_Lenth - 1;  
    
    for(Count = 0; Count < b_Lenth ;Count++)
    {    	
            Output_Data += (*(b + Count)) *
                            (*(Input_Data - Count));                         
    }         
    
    return (double)Output_Data;
}





int main(void)
{
    double w_p = pi/10;
    double w_s = pi/5;
    double w_c = (w_s + w_p)/2;
    printf("w_c =  %f \n" , w_c);
    
    int N = 0;  
    N = (int) ((6.6*pi)/(w_s - w_p) + 0.5);
    if(0 == N%2) N++;    
    printf("N =  %d \n" , N);    
  
    double *Impulse_Response;        
    Impulse_Response = (double *) malloc(sizeof(double)*N);  
    memset(Impulse_Response,
           0,
           sizeof(double)*N);
           
    Unit_Impulse_Response(N,w_c,
                          Hamming,
                          Impulse_Response);       

    double *Input_Buffer;
    double Output_Data = 0;
    Input_Buffer = (double *) malloc(sizeof(double)*N);  
    memset(Input_Buffer,
           0,
           sizeof(double)*N);
    int Count = 0;
    
    FILE *fs; 
    fs=fopen("FIR_Data.txt","w"); 
    
    while(1)
    {   
    	if(Input_Data[Count] == 55) break;
    	 
    	Save_Input_Date (Input_Data[Count],
        	         N,
                	 Input_Buffer);
       
        Output_Data = Real_Time_FIR_Filter(Impulse_Response,
                                           N,
                                           Input_Buffer);
                 		
          
        fprintf(fs,"%lf,",Output_Data);
        //if(((Count+1)%5 == 0)&&(Count != 0))  fprintf(fs,"\r\n"); 
   
        Count++;
    }
          
    /*---------------------display--------------------------------     
    for(Count = 0; Count < N;Count++)
    {
        printf("%d %lf  ",Count,*(Input_Buffer+Count));
        if(((Count+1)%5 == 0)&&(Count != 0)) printf("\n");        	 
    }      
    */  
    
    fclose(fs);
    printf("Finish \n");
    return (int)0;
}
           

3.頻響的問題

        按照上面程式,參數如下設定。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

        運作程式,我們就實作了一個FIR濾波器。我們使用Matlab做出其頻響。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

        好了,這裡可以看出,從其幅度特性看來,我們确實實作了一個低通濾波器。但是,相位特性就比較奇怪(為了友善看出問題,我已經進行了解卷繞,至于什麼是解卷繞,為什麼要解卷繞,之後會說)。

       那麼,問題來了!按照道理來說,這個FIR濾波器應該是擁有線性相位特性的,但是為什麼這裡的線性相位特性确不是一條直線!在接近于阻帶起始頻率的地方,有什麼會震蕩?這個問題,之後再解決吧。

      [數字信号處理]相位特性解卷繞   <-----------戳我

4.實際的濾波效果

       為了驗證效果,我們輸入實際的信号看看。這裡,我們選擇采樣周期為10ms(0.1ms,2014.10.3日修正),那麼,其通帶頻率和阻帶起始頻率為

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

       為了驗證其性質,我選擇了1KHz和3KHz的頻率混合,最終輸出。輸入的波形如下。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

      其輸出波形是如下。

[數字信号處理]機關沖擊響應與頻響以及FIR實作代碼(C語言)

        紅色的“+”是Matlab計算的結果,黑的o是我用C語言實作的濾波器的計算結果,藍的*号是1KHz的信号,也就是希望的輸出。可以看出,這個濾波器有一定的延遲,但是濾波效果還是不錯。

       部落格位址:http://blog.csdn.net/thnh169/ 

繼續閱讀