DSP

低通滤波器差分方程用C语言实现的问题

2019-07-15 18:56发布

站内问答 / DSP
14656 6 1598
本帖最后由 hanl 于 2011-7-21 10:07 编辑


我用C语言实现一个butterworth低通滤波器，现在遇到问题，please help me！
实现过程：
1、采样频率10KHz，通带边界频率2KHz，通带波纹小于1dB，阻带边界频率3KHz，阻带衰减大于40dB。
2、用matlab求出滤波器系数b，a
  1. >> Fp = 2000;
  2. >> Fs = 3000;
  3. >> F = 10000;
  4. >> T = 1/F;
  5. >> Ap = 1;
  6. >> As = 40;
  7. >> Wp = Fp*T*2;
  8. >> Ws = Fs*T*2;
  9. >> [n,Wc] = buttord(Wp,Ws,Ap,As);
  10. >> [b,a] = butter(n,Wc)

  12. b =

  14. 0.0021 0.0186 0.0745 0.1739 0.2609 0.2609 0.1739 0.0745 0.0186 0.0021


  17. a =

  19. 1.0000 -1.0893 1.6925 -1.0804 0.7329 -0.2722 0.0916 -0.0174 0.0024 -0.0001

  21. >>
复制代码3、差分方程： y(0) = b(0)*x(0) + b(1)*x(1) +...+ b(9)*x(9) - {a(1)*y(1) + a(2)*y(2) +...+ a(9)*y(9)}
4、程序：
  1. double a[L] = {1.0000, -1.0893, 1.6925, -1.0804, 0.7329, -0.2722, 0.0916, -0.0174, 0.0024, -0.0001};
  2. double b[L] = {0.0021, 0.0186, 0.0745, 0.1739, 0.2609, 0.2609, 0.1739, 0.0745, 0.0186, 0.0021};

  4. read(&input); //输入范围：-10V～+10V ，量化范围 0～65535
  5. x[0] = input-32767;  
  6. for (i=0; i<=L-1; i++)
  7. {
  8. b_sum = b[i]*x[i] + b_sum;
  9. a_sum = a[i+1]*y[i+1] + a_sum;
  10. }
  11. y[0] = b_sum - a_sum;
  12. output = y[0];
  13. for (i=L-1; i>=1; i--)
  14. {
  15. x[i] = x[i-1];
  16. y[i] = y[i-1];
  17. }
  18. write(output);
复制代码5、测试结果
   a、无输入
Screenshot1.png   （单位：0.1ms）

b、输入1KHz正弦波，峰峰值 4V
Screenshot2.png   （单位：0.1ms）

到此，不知道哪里出了问题，请高人赐教
友情提示: 此问题已得到解决，问题已经关闭,关闭后问题禁止继续编辑，回答。
该问题目前已经被作者或者管理员关闭, 无法添加新回复
6条回答
火奴小子
1楼-- · 2019-07-16 00:45
传递函数应该是AY=BX吧，所以你的差分方程也好像不对吧！！！！！！！a，b系数反了吧
hanl
2楼-- · 2019-07-16 02:39
 精彩回答 2  元偷偷看……
fylwzh
3楼-- · 2019-07-16 06:02
{:1:}{:1:}{:1:}
liugangcumt
4楼-- · 2019-07-16 11:18
你好，问题找到吗？在学习数字滤波。请指教。
jindongli
5楼-- · 2019-07-16 12:36
差分方程是Z域的，你要求的是时域的。你如果要用差分方程就要Z逆变换，要不就用卷积。
jindongli
6楼-- · 2019-07-16 16:14
滤波函数                          0.15
                        H(z) = ----------------------
                                        1 - 0.8*z^-1
1.        b = 0.15 ;  
2.        a = [1 -0.8];  
3.          
4.        n = [0:100];  
5.        x = 2*sin(0.05*pi*n) + 2*randn(1, 101); %滤波前序列  
6.          
7.        imp = [1; zeros(100, 1)];  
8.        h = filter(b, a, imp); % h为冲激函数的响应（时域）  
9.        w = conv(x, h); % 做卷积  
10.        yc = w(1:101);  %滤波后序列   
11.        z= filter(b, a, x); % 进行滤波   
12.        xlabel('n');  
13.        ylabel('x y yc');  
14.        plot(n, x,'r',n,z ,'b', n,yc,'m');  
15.        grid;  

思想很简单，首先是用冲击函数 Imp = [1 ; zeros(100,1)] 取求出系统的冲激响应，然后再与输入进行卷积
运算结果：
  

一周热门 更多>