[数字信号处理]单位冲击响应与频响以及FIR实现代码（C语言）

1.单位冲击响应与频响

        就如同之前所说的一样，使用下图所示的单位冲击响应，所设计的滤波器，是无法实现的。

         现在，让我们看看其这个滤波器的频响。所谓频响，就是计算其单位冲击响应的离散时间傅里叶变换，

       我们可以看出，这个滤波器的频响的计算结果是实数，并没有虚数部分。也就是，其相位谱一直是0，也就意味着，这个滤波器输入与输出之间没有延迟，这种相位特性称为0延迟相位特性。

        但是，这个滤波器无法是无法实现的。我们实际计算一下该滤波器的输入输出，就可以明白。

        这个滤波器在计算的过程中，需要过去的值和未来的值。未来的值是不可预测的，所以，这个滤波器无法实现。为了使得这个滤波器可以实现，我们只好移动其单位冲击响应，使得其不再需要未来的值。比如，就像下面这样移动。

        这样的话，这个滤波器就可以实现了，我们只需要记录其40个过去的输入值和现在的输入值。但是，由于移动了其单位冲击响应，会不会对频响产生什么影响呢，我们来看看。

       为了更好的说明问题，L去代替之前例子中的20。

       移动之后频响，我们根据上面式子可以得出两个结论：1，移动不会对幅度谱产生影响。2,，移动会对相位产生一个延迟，这个延迟主要取决于移动的长度，移动的长度越长，延迟越大。但是，这个移动是线性的。

       因此，我们把这个移动的相位特性称为，线性相位特性。到这里，我们移动后的，因果的，可实现的滤波器的单位冲击响应，如下所示。

2.窗函数实现的FIR滤波器代码（C语言）

[cpp]  view plain copy

1. #include <stdio.h>  
2. #include <math.h>  
3. #include <malloc.h>  
4. #include <string.h>  
5.   
6.   
7. #define   pi         (3.1415926)  
8.   
9. /*-------------Win Type----------------*/  
10. #define   Hamming    (1)  
11.   
12.   
13.   
14. double Input_Data[] =   
15. {  
16. 0.000000 , 0.896802 , 1.538842 , 1.760074 ,  1.538842 ,  1.000000 ,  0.363271 , -0.142040 , -0.363271 , -0.278768,  
17. 0.000000 , 0.278768 , 0.363271 , 0.142020 , -0.363271 , -1.000000 , -1.538842 , -1.760074 , -1.538842 , -0.896802,  
18. 0.000000 , 0.896802 , 1.538842 , 1.760074 ,  1.538842 ,  1.000000 ,  0.363271 , -0.142040 , -0.363271 , -0.278768,  
19. 0.000000 , 0.278768 , 0.363271 , 0.142020 , -0.363271 , -1.000000 , -1.538842 , -1.760074 , -1.538842 , -0.896802,  
20. 0.000000 , 0.896802 , 1.538842 , 1.760074 ,  1.538842 ,  1.000000 ,  0.363271 , -0.142040 , -0.363271 , -0.278768,  
21. 0.000000 , 0.278768 , 0.363271 , 0.142020 , -0.363271 , -1.000000 , -1.538842 , -1.760074 , -1.538842 , -0.896802,  
22. 0.000000 , 0.896802 , 1.538842 , 1.760074 ,  1.538842 ,  1.000000 ,  0.363271 , -0.142040 , -0.363271 , -0.278768,  
23. 0.000000 , 0.278768 , 0.363271 , 0.142020 , -0.363271 , -1.000000 , -1.538842 , -1.760074 , -1.538842 , -0.896802,  
24. 0.000000 , 55  
25. };  
26.   
27.   
28.   
29.   
30. double sinc(double n)  
31. {  
32.     if(0==n) return (double)1;  
33.     else return (double)(sin(pi*n)/(pi*n));  
34. }  
35.    
36. int Unit_Impulse_Response(int N,double w_c,  
37.                           int Win_Type,  
38.                           double *Output_Data)  
39. {  
40.     signed int Count = 0;   
41.     
42.     for(Count = -(N-1)/2;Count <= (N-1)/2;Count++)  
43.     {  
44.         *(Output_Data+Count+((N-1)/2)) = (w_c/pi)*sinc((w_c/pi)*(double)(Count));  
45.         //printf("%d %lf  ",Count+((N-1)/2)+1,*(Output_Data+Count+((N-1)/2)));  
46.         //if(Count%4 == 0) printf("\n");  
47.     }     
48.       
49.       
50.     switch (Win_Type)  
51.     {  
52.           
53.         case Hamming:   printf("Hamming \n");  
54.                         for(Count = -(N-1)/2;Count <= (N-1)/2;Count++)  
55.                 {  
56.                     *(Output_Data+Count+((N-1)/2)) *= (0.54 + 0.46 * cos((2*pi*Count)/(N-1)));  
57.                     //printf("%d %lf  ",Count+((N-1)/2)+1,*(Output_Data+Count+((N-1)/2)));  
58.                     //if(((Count+1)%5 == 0)&&(Count != -20)) printf("\n");  
59.             }   
60.                         break;  
61.                           
62.                           
63.         default:        printf("default Hamming \n");  
64.                         for(Count = -(N-1)/2;Count <= (N-1)/2;Count++)  
65.                 {  
66.                     *(Output_Data+Count+((N-1)/2)) *= (0.54 + 0.46 * cos((2*pi*Count)/(N-1)));  
67.                     //printf("%d %lf  ",Count+((N-1)/2)+1,*(Output_Data+Count+((N-1)/2)));  
68.                     //if(((Count+1)%5 == 0)&&(Count != -20)) printf("\n");  
69.             }   
70.           
71.                         break;  
72.     }  
73.       
74.     return (int)1;  
75. }  
76.   
77.   
78. void Save_Input_Date (double Scand,  
79.                       int    Depth,  
80.                       double *Input_Data)  
81. {  
82.     int Count;  
83.     
84.     for(Count = 0 ; Count < Depth-1 ; Count++)  
85.     {  
86.         *(Input_Data + Count) = *(Input_Data + Count + 1);  
87.     }  
88.       
89.     *(Input_Data + Depth-1) = Scand;  
90. }  
91.   
92.   
93.   
94. double Real_Time_FIR_Filter(double *b,  
95.                             int     b_Lenth,  
96.                             double *Input_Data)  
97. {      
98.     int Count;  
99.     double Output_Data = 0;  
100.       
101.     Input_Data += b_Lenth - 1;    
102.       
103.     for(Count = 0; Count < b_Lenth ;Count++)  
104.     {         
105.             Output_Data += (*(b + Count)) *  
106.                             (*(Input_Data - Count));                           
107.     }           
108.       
109.     return (double)Output_Data;  
110. }  
111.   
112.   
113.   
114.   
115.   
116. int main(void)  
117. {  
118.     double w_p = pi/10;  
119.     double w_s = pi/5;  
120.     double w_c = (w_s + w_p)/2;  
121.     printf("w_c =  %f \n" , w_c);  
122.       
123.     int N = 0;    
124.     N = (int) ((6.6*pi)/(w_s - w_p) + 0.5);  
125.     if(0 == N%2) N++;      
126.     printf("N =  %d \n" , N);      
127.     
128.     double *Impulse_Response;          
129.     Impulse_Response = (double *) malloc(sizeof(double)*N);    
130.     memset(Impulse_Response,  
131.            0,  
132.            sizeof(double)*N);  
133.              
134.     Unit_Impulse_Response(N,w_c,  
135.                           Hamming,  
136.                           Impulse_Response);         
137.   
138.     double *Input_Buffer;  
139.     double Output_Data = 0;  
140.     Input_Buffer = (double *) malloc(sizeof(double)*N);    
141.     memset(Input_Buffer,  
142.            0,  
143.            sizeof(double)*N);  
144.     int Count = 0;  
145.       
146.     FILE *fs;   
147.     fs=fopen("FIR_Data.txt","w");   
148.       
149.     while(1)  
150.     {     
151.         if(Input_Data[Count] == 55) break;  
152.            
153.         Save_Input_Date (Input_Data[Count],  
154.                      N,  
155.                      Input_Buffer);  
156.          
157.         Output_Data = Real_Time_FIR_Filter(Impulse_Response,  
158.                                            N,  
159.                                            Input_Buffer);  
160.                           
161.             
162.         fprintf(fs,"%lf,",Output_Data);  
163.         //if(((Count+1)%5 == 0)&&(Count != 0))  fprintf(fs,"\r\n");   
164.      
165.         Count++;  
166.     }  
167.             
168.     /*---------------------display--------------------------------      
169.     for(Count = 0; Count < N;Count++) 
170.     { 
171.         printf("%d %lf  ",Count,*(Input_Buffer+Count)); 
172.         if(((Count+1)%5 == 0)&&(Count != 0)) printf("\n");            
173.     }       
174.     */    
175.       
176.     fclose(fs);  
177.     printf("Finish \n");  
178.     return (int)0;  
179. }  

3.频响的问题

        按照上面程序，参数如下设定。

   

        运行程序，我们就实现了一个FIR滤波器。我们使用Matlab做出其频响。

        

        好了，这里可以看出，从其幅度特性看来，我们确实实现了一个低通滤波器。但是，相位特性就比较奇怪（为了方便看出问题，我已经进行了解卷绕，至于什么是解卷绕，为什么要解卷绕，之后会说）。

       那么，问题来了！按照道理来说，这个FIR滤波器应该是拥有线性相位特性的，但是为什么这里的线性相位特性确不是一条直线！在接近于阻带起始频率的地方，有什么会震荡?这个问题，之后再解决吧。

      [数字信号处理]相位特性解卷绕   <-----------戳我

4.实际的滤波效果

       为了验证效果，我们输入实际的信号看看。这里，我们选择采样周期为10ms，那么，其通带频率和阻带起始频率为

       为了验证其性质，我选择了1KHz和3KHz的频率混合，最终输出。输入的波形如下。

      其输出波形是如下。

        红色的“+”是Matlab计算的结果，黑的o是我用C语言实现的滤波器的计算结果，蓝的*号是1KHz的信号，也就是希望的输出。可以看出，这个滤波器有一定的延迟，但是滤波效果还是不错。