原创数字滤波器matlab辅助设计

 2022-6-5 07:38  2671 18 11 分类: MCU/ 嵌入式文集: matlab

仍然使用FDATOOL设计滤波器，当前设计一个数字带通滤波器。至于用的是冲击响应不变法，还是其它的方法。暂时不考虑。FIR 需要的阶数太多，也不考虑。使用IIR滤波，线性相位就不要想了。可以选巴特沃兹（最大平整度），或切比雪夫（最大陡降特性。）发现在相同的性能下切比雪夫需要的阶数少。

生成的传递函数是按照多个二阶单元级联。系统提供 SOS（Second Order Section）也可以称为“救命”矩阵。其思路是将高阶传递函数分解为多个稳定的二阶函数级联，保持系统稳定。因为使用的是IIR （无限冲击响应）构成，注定其是非线性相位。会导致滤波后的信号波形畸变。导出滤波器参数文件.

每个SOS参数可以使用直接II型，实现如下：

编写一段代码测试一下。

clear
close all
Fs = 30000; % Sampling Frequency
% 20-500 30000sps butter 4 order
sec = 2;
SOS = [1 0 -1 1 -1.9941159332576488 0.99413480892694572;1 0 -1 1 -1.8630805194585167 0.87259004111123051];
G = [0.048551094455762382 0.048551094455762382 1];
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 = [1 0 -1 1 -1.9941159332576488 0.99413480892694572;1 0 -1 1 -1.8630805194585167 0.87259004111123051];
G = [0.048551094455762382 0.048551094455762382 1];
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[SEC][6] = {
{1, 0, -1, 1, -1.9528158412163148, 0.95399450820943588},
{1, 0, -1, 1, -0.92310027138083539, 0.34076352170725532}};
const float G1[SEC+1] = {0.31856548677578983,0.31856548677578983,1};
float BandPassFilter(const float SOS[SEC][6],const float G[SEC+1],float x){
static float w[SEC][6];
uint8_t k;
static float ym[SEC+1];
ym[0] = x;
for(k = 0;k < SEC;k++){
w[k][0] = ym[k] * G[k] - w[k][1]*SOS[k][4] - w[k][2]*SOS[k][5];
ym[k+1] = w[k][0]*SOS[k][0]+w[k][1]*SOS[k][1]+w[k][2]*SOS[k][2];
w[k][2] = w[k][1];
w[k][1] = w[k][0];
}
return ym[SEC]*G[SEC];
}