标签: SOS

仍然使用 FDATOOL 设计滤波器，当前设计一个数字带通滤波器。至于用的是冲击响应不变法，还是其它的方法。暂时不考虑。 FIR 需要的阶数太多，也不考虑。使用 IIR 滤波，线性相位就不要想了。可以选巴特沃兹（最大平整度），或切比雪夫（最大陡降特性。）发现在相同的性能下切比雪夫需要的阶数少。 生成的传递函数是按照多个二阶单元级联。系统提供 SOS （ Second Order Section ）也可以称为“救命”矩阵。其思路是将高阶传递函数分解为多个稳定的二阶函数级联，保持系统稳定。因为使用的是 IIR （无限冲击响应）构成，注定其是非线性相位。会导致滤波后的信号波形畸变。 导出滤波器参数文件 . % Generated by MATLAB(R) 8.4 and the Signal Processing Toolbox 6.22. % Generated on: 04-Jun-2022 20:49:14 % Coefficient Format: Decimal % Discrete-Time IIR Filter (real) % ------------------------------- % Filter Structure : Direct-Form II, Second-Order Sections % Number of Sections : 2 % Stable : Yes % Linear Phase : No SOS Matrix: 1 0 -1 1 -1.9528158412163148 0.95399450820943588 1 0 -1 1 -0.92310027138083539 0.34076352170725532 Scale Values: 0.31856548677578983 0.31856548677578983 每个 SOS 参数可以使用直接 II 型，实现如下： 编写一段代码测试一下。 clear close all Fs = 30000; % Sampling Frequency % 20-500 30000sps butter 4 order sec = 2; SOS = ; G = ; len = Fs/10; % Êý¾Ý³¤¶ÈΪ 100mS t = 1:len; n = 300*sin(2*pi*15*t/Fs + pi/6); % 50hz noise s = 3000*sin(2*pi*200*t/Fs + pi/5); % original square wave signal ng = 500*randn(1,len); x = s + n + + ng + 1000; fftx = abs(fft(x,len)); f = 0:Fs/len:(len-1)*(Fs/len); figure(1) plot(f(1:len/20),fftx(1:len/20)) title('sinc fuction spectrum') w = zeros(sec,3); ym = zeros(1,sec+1); for i = 1:len ym(1) = x(i); for k = 1:sec w(k,1) = ym(k) * G(k) - w(k,2) * SOS(k,5) - w(k,3) * SOS(k,6) ym(k+1) = w(k,1)*SOS(k,1)+w(k,2)*SOS(k,2)+w(k,3)*SOS(k,3); w(k,3) = w(k,2); w(k,2) = w(k,1); end y(i) = ym(sec+1)*G(sec+1); end figure(2) plot(t,s,'r',t,x,'c',t,y,'.b') legend('signal','org','bandpass'); 待滤波信号的频谱特性 可以看到 org 信号上有噪声，有直流偏置和低频交流干扰。经过带通滤波后，去除了直流，去除了噪声，基本还原原始信号 signal 。刚开始时并不稳定需要等一段时间才能达到稳定。 可见，此段代码可以将一个输入序列 x 转换为一个输出序列 y ，完成对 x 序列的滤波。 按照生成的结构图，编写 m 文件的解释： sec = 2; % 定义级数 SOS = ; G = ; w = zeros(sec,3); %中间w 变量 ym = zeros(1,sec+1); %每一级的输出 for i = 1:len % 对每一个输入 x ym(1) = x(i); for k = 1:sec % 分极计算 w(k,1) = ym(k) * G(k) - w(k,2) * SOS(k,5) - w(k,3) * SOS(k,6) ym(k+1) = w(k,1)*SOS(k,1)+w(k,2)*SOS(k,2)+w(k,3)*SOS(k,3); w(k,3) = w(k,2); w(k,2) = w(k,1); end y(i) = ym(sec+1)*G(sec+1); %得到输出 end 翻译成C语言 // 20-500hz butterworth 4 order 2 section 3750sps #define SEC 2 const float SOS1 = { {1, 0, -1, 1, -1.9528158412163148, 0.95399450820943588}, {1, 0, -1, 1, -0.92310027138083539, 0.34076352170725532}}; const float G1 = {0.31856548677578983,0.31856548677578983,1}; float BandPassFilter(const float SOS ,const float G ,float x){ static float w ; uint8_t k; static float ym ; ym = x; for(k = 0;k < SEC;k++){ w = ym * G - w *SOS - w *SOS ; ym = w *SOS +w *SOS +w *SOS ; w = w ; w = w ; } return ym *G ; } 每采集到一个数据时，以输入数据作为参数，调用 BandPassFilter （），得到一个滤波输出。相当于实时处理。 经过滤波后的信号 经过滤波 + 陷波后的信号