【滤波器】matlab仿真巴特沃斯滤波器并在stm32进行运行

进击的wayne

已于 2024-05-06 22:17:40 修改

阅读量840

点赞数 20

文章标签： matlab stm32 c语言

于 2024-05-06 21:52:37 首次发布

本文链接：https://blog.csdn.net/qq_34975519/article/details/138509444

版权

【滤波器】matlab仿真巴特沃斯滤波器并在stm32进行运行

摘要

本文记录学习数字滤波器的过程，大致包含五个部分：数字滤波器的理解，fdatool的使用，窗函数，matlab生成数字滤波器，通过stm32f103对滤波器进行仿真。其中对fdatool的使用和stm32的使用在其他文章中进行了简述。末文附上学习过程中的参考链接，非常感谢各位大佬的学习分享。

1 滤波器的分类

1.1 数字滤波器和模拟滤波器

区别：

（1）模拟滤波器处理的是连续时间信号，而数字滤波器处理的是离散时间信号。
（2）模拟滤波器使用的是模拟电路技术，其输入信号和输出信号都是模拟信号。数字滤波器使用的是数字信号处理技术，其输入信号和输出信号都是数字信号。
（3）数字滤波器关于Fs/2频率是翻转的，也就是对称的，而模拟滤波器则不具备这样的特性

实现方式：

数字滤波器通过算法在数字处理器上实现。
模拟滤波器通过电子元件（如电容、电阻、电感）在模拟电路中实现。

应用范围：

数字滤波器广泛用于数字信号处理、通信系统、音频处理等领域。 模拟滤波器通常用于模拟电路、音频放大器、射频电路等领域。

数字滤波器：

数字滤波器有频域滤波器和时域滤波器之分，常见的时域滤波器又分为有限冲激响应FIR滤波器和无限冲激响应IIR滤波器。FIR滤波器的具有线性相位特性而IIR滤波器是非线性相位。
数字滤波器可以分为低通滤波器、带通滤波器、高通滤波器、带阻滤波器、椭圆滤波器、抽头滤波器、滤波器组等。
任何一个系统都可以被认为是一种滤波器。滤波器又具备着“幅频相应”和“相频相应”的概念，所谓的幅频相应就是对输入信号的某个频率的分量进行了增强，而相频相应是对某个频率分量信号的时延的。`<img src="image/filter/1714546583507.png" alt="Figure 1" />`

1.2 策略不同的滤波器

(1) 常见的IIR滤波器包括一下几种，

滤波器名称	特性
巴特沃斯滤波器	是在通频带内实现尽可能平坦的频率响应，即在通带内没有波纹。这意味着它具有最大平坦度的通频带响应。
切比雪夫滤波器	给定通频带内实现最小的通带波纹，即在通带内具有最小的振荡，同时在通频带内具有更快的幅度下降速度，但通带内会有一定程度的波纹，这是为了在通带内获得更快的过渡。
贝塞尔滤波器	现最小相位失真，即在时域上尽可能保持信号的形状，避免引入额外的相位延迟，有相对较慢的幅度下降速度，但在时域上具有最小的相位失真，因此在一些对信号相位要求较高的应用中很有用
椭圆滤波器

在这里插入图片描述

IIR 表达公式

$H(z)=\frac{\sum_{i=1}^{N}b_iz^-i}{1-\sum_{l=1}^{N}a_lz^-l} \\ y[n]=\sum_{i=1}^{N}x(n-1)b(i)+\sum_{l=1}^{N}y(n-l)a(l) \\ 2Biquad直接I型：y[n]=b_1*x[n]+b_2x[n-1]+b_3x[n-2]+a_1*y[n-1]+a_2*y[n-2] \\ 2Biquad直接II型:y[n]_1=X[n]-a_1*X[n-1]-a_1*X[n-2] \\ y[n]=b_1*y[n]_1+b_2*y[n-1]_1+b_3*y[n-2]_1$

系统框图

在这里插入图片描述

(2) Fir滤波器学习链接：
参考学习链接：https://blog.csdn.net/yishuihanq/article/details/131711885

1.3高通低通滤波器设计

学习链接：https://blog.csdn.net/ccsss22/article/details/136788020

2 窗函数

信号加窗是为了避免频谱泄露，频谱泄漏是指在进行离散傅里叶变换时,由于信号不是周期性的,导致频域上的能量泄漏到相邻频率上的现象。这会导致变换结果失真,影响信号分析和处理的准确性。

学习链接：https://blog.csdn.net/weixin_63595975/article/details/134277281

3 fdatool使用

另一篇文章介绍：https://mp.csdn.net/mp_blog/creation/success/138510176

4 matlab 生成butterworth低通滤波器实验

实验：模拟生成1~10hz的正选信号，并分别通过fdatool生成和matlab函数生成的方式生成1hz的巴特沃斯滤波器。对最终结果进行显示。

代码解析：

close all;

fs=1024;% 采样频率
t=0:1/fs:6-1/fs; % 横坐标 采样点数=(6-1/fs)
freq_Matrix=[1;2;3;4;5;6;7;8;9;10];

%生成原始信号
y=zeros(size(t));
for i=1:length(freq_Matrix)
    y=y+sin(2*pi*freq_Matrix(i)*t);
end

%对信号进行fft
f=(0:length(y)-1)*(fs/length(y));
fft_y=fft(y);

% 通过butter函数生成butterworth滤波器
fc_low=1;
wn=2*fc_low/fs;
[b,a]=butter(8,wn,'low');
freqz(b,a,[],fs);% 生成幅频相位曲线

% 使用代码生成的巴特沃斯低通器进行滤波
y_low=filter(b,a,y);
y_1hz=sin(2*pi*t)+sin(4*pi*t);

% 通过fdatool生成的巴特沃斯低通器进行滤波
y_low_fda=filter(bwfilterlowpass1hz,y);

%figure显示
figure('Name','原始信号','NumberTitle','on');
subplot(3,2,1);plot(t,y/length(freq_Matrix));
subplot(3,2,2);plot(f(1:length(y)/50),abs(fft_y(1:length(y)/50)));
subplot(3,2,3);plot(t,y_low);
subplot(3,2,4);plot(t,y_1hz);
subplot(3,2,5);plot(t,y_low_fda);
subplot(3,2,6);plot(t,y_1hz);% 结果相同

实验现象：

6 stm32f103c8t6实验

（1）实验失败案例

实验：通过stm32f103c模拟生成1~10hz的正选信号，并对信号进行butterworth低通滤波，通过keil的逻辑分析仪对数据进行回显分析。

参数设置：

直接I型参数：
SOS Matrix:
1  2  1  1  -1.9952772372946446  0.99531479853184868
1  2  1  1  -1.9886888350211553  0.98872627223121401

Scale Values:
0.0000093903093010300073
0.0000093593025146089099

直接II型参数：
SOS Matrix:
1  2  1  1  -1.9952772372946446  0.99531479853184868
1  2  1  1  -1.9886888350211553  0.98872627223121401

Scale Values:
0.0000093903093010300073
0.0000093593025146089099

代码解析：

float g_lowpass_1hz_IIRCoeffs_I[6 * 2] = {1, 2, 1, 1, -1.9952772372946446, 0.99531479853184868,
                                          1, 2, 1, 1, -1.9886888350211553, 0.98872627223121401};

float g_lowpass_1hz_Scalevalue_I[2] = {0.0000093903093010300073, 0.0000093593025146089099};

float g_lowpass_1hz_IIRCoeffs_II[6 * 2] = {1, 2, 1, 1, -1.9952772372946446, 0.99531479853184868,
                                           1, 2, 1, 1, -1.9886888350211553, 0.98872627223121401};

float g_lowpass_1hz_Scalevalue_II[2] = {0.093903093010300073, 0.093593025146089099};

int16_t sin_1HZ, sin_2HZ, sin_3HZ, sin_4HZ, sin_5HZ, sin_6HZ;
float sin_all;
float pDelayI[8] = {0}; // Xn1 Xn2 Yn1 Yn2 二阶的话就有两个 2*4
float pDelayII[6] = {0};
volatile float _wave_1hz = 0; // 1hz的原始信号
volatile float _waveI, _waveII = 0;
volatile float _wave_1hz_frI = 0; // 直接I型 filter result信号
float _wave_1hz_frII = 0;         // 直接II型 filter result信号

/**
 *  @brief butterworth 滤波器直接一型计算公式
 *          y[n]=b1*x[n]+b2*x[n-1]+b3*x[n-2]+a1*y[n-1]+a2*y[n-2]
 *          y[n]=scale*(b1*x[n]+b2*x[n-1]+b3*x[n-2])-a1*y[n-1]-a2*y[n-2]
 *  @param  p_IIRCoeffs IIR滤波器系数
 *  @param  p_Scalevalue 缩放系数
 *  @param  p_Delay 保存的状态变量 直接1型四个时延单位
 *  @param  _Input 此刻输入值
 *  @param  _stage // 有几个biquad 一个biquad对应一个session
 *  参考链接：https://blog.csdn.net/ether_7/article/details/119795721
 *  @return
 */
float butterWorth2Order_I_Float(float *p_IIRCoeffs, float *p_Scalevalue, float *p_Delay, float _input, int16_t _stage)
{
    float output;
    int16_t _coefStep, _delayStep;
    float temp[2] = {0};
    float _scale = 0;
    for (int16_t i = 0; i < _stage - 1; i++)
    {
        _coefStep = i * 6;
        _delayStep = i * 4;
        _scale = p_Scalevalue[i];
        temp[i] = _scale * (p_IIRCoeffs[_coefStep] * _input + p_IIRCoeffs[_coefStep + 1] * p_Delay[_delayStep] +
                            p_IIRCoeffs[_coefStep + 2] * p_Delay[_delayStep + 1]);

        output = temp[i] - p_IIRCoeffs[_coefStep + 4] * p_Delay[_delayStep + 2] -
                 p_IIRCoeffs[_coefStep + 5] * p_Delay[_delayStep + 3];

        p_Delay[_delayStep] = _input;                      // Xn
        p_Delay[_delayStep + 1] = p_Delay[_delayStep];     // X[n-1]=X[n]
        p_Delay[_delayStep + 2] = output;                  // Yn
        p_Delay[_delayStep + 3] = p_Delay[_delayStep + 2]; // Y[n-1]=Y[n]
        _input = output;                                   // 下一级的输入等于上一级的输出
    }
    return output;
}
/**
 *  @brief butterworth 滤波器直接II型计算公式
 *          y[n]=scale*(b1*x[n]+b2*x[n-1]+b3*x[n-2])-a1*y[n-1]-a2*y[n-2]
 *  @param 与直接I型相同 直接II型3个时延单位
 *  @return
 */
volatile int16_t tes = 0;
float butterWorth2Order_II_Float(float *p_IIRCoeffs, float *p_Scalevalue, float *p_Delay, float _input, uint16_t _stage)
{
    float output = 0;
    uint16_t _coefStep, _delayStep;
    float _scale = 0;
    float b0, b1, b2;
    float a1, a2;
    for (uint16_t i = 0; i < _stage - 1; i++)
    {
        _coefStep = i * 6;
        _delayStep = i * 3;
        _scale = p_Scalevalue[i];
        p_Delay[_delayStep] = _scale * _input - p_IIRCoeffs[_coefStep + 4] * p_Delay[_delayStep + 1] -
                              p_IIRCoeffs[_coefStep + 5] * p_Delay[_delayStep + 2];

        output = p_IIRCoeffs[_coefStep] * p_Delay[_delayStep] + p_IIRCoeffs[_coefStep + 1] * p_Delay[_delayStep + 1] +
                 p_IIRCoeffs[_coefStep + 2] * p_Delay[_delayStep + 2];

        p_Delay[_delayStep + 2] = p_Delay[_delayStep + 1]; // X[n-2]=X[n-1]
        p_Delay[_delayStep + 1] = p_Delay[_delayStep];     // X[n-1]=X[N]

        _input = output; // 下一级的输入等于上一级的输出
    }

    return output;
}
/**
 *  @brief butterworth 1hz 滤波器 c代码实现
 *  @param
 *  @return
 */
void Exp_BwFilter_Lowpass_1HZ_Float(float *result)
{

    uint32_t divide = 0;

    while (1)
    {
        // 1step: 生成1~6hz的信号
        sin_1HZ = Generate_sin(divide, 1);
        sin_2HZ = Generate_sin(divide, 2);
        sin_3HZ = Generate_sin(divide, 3);
        sin_4HZ = Generate_sin(divide, 4);
        sin_5HZ = Generate_sin(divide, 5);
        sin_6HZ = Generate_sin(divide, 6);
        sin_all = sin_1HZ + sin_2HZ + sin_3HZ + sin_4HZ + sin_5HZ + sin_6HZ;
        _waveI = Generate_Wave_ByAmP(100, sin_all);
        _waveII = _waveI;
        _wave_1hz = Generate_Wave_ByAmP(100, sin_1HZ);

        divide++;

        // // 2step: 进行巴特沃斯滤波
        _wave_1hz_frI =
            butterWorth2Order_I_Float(g_lowpass_1hz_IIRCoeffs_I, g_lowpass_1hz_Scalevalue_I, pDelayI, _waveI, 2);
        _wave_1hz_frII =
            butterWorth2Order_II_Float(g_lowpass_1hz_IIRCoeffs_II, g_lowpass_1hz_Scalevalue_II, pDelayII, _waveII, 2);
    }
}

实验现象：

一开始的时候这个实验失败了，不是代码写的有问题，是因为keil并不直接支持显示double类型数据。double类型数据通常用于表示浮点数，而逻辑分析仪更适用于观察和分析数字信号的逻辑状态，比如高电平和低电平。如果你想观察浮点数的值，你可以通过在代码中使用printf函数将double类型数据转换为字符串，然后将字符串输出到串口或其他可以显示文本的终端设备，如串口调试助手。然后，你可以通过逻辑分析仪观察串口输出的数据信号。
所以为了保证现象符合理论，将实验过程中的数据以牺牲精度为代价都改变采用float类型。然后用keil中进行debug 将变量加到analyzer进行查看
此外，感觉这个直接I型的结果有问题，但是调了好久也还是这样，如果有发现问题的大哥们能帮我提出来一下。

在这里插入图片描述

参考链接：

markdown公式大全： https://blog.csdn.net/jzj_c_love/article/details/122279703
fir滤波器学习：https://blog.csdn.net/yishuihanq/article/details/131711885
窗函数：https://blog.csdn.net/weixin_63595975/article/details/134277281
高通低通滤波器设计：https://blog.csdn.net/ccsss22/article/details/136788020
matlab函数查询网：https://ww2.mathworks.cn/
巴特沃斯滤波器c实现代码：https://blog.csdn.net/ether_7/article/details/119795721

进击的wayne

关注

20
点赞
踩
11

收藏

觉得还不错? 一键收藏
0
评论
【滤波器】matlab仿真巴特沃斯滤波器并在stm32进行运行

本文记录学习数字滤波器的过程，大致包含五个部分：数字滤波器的理解，fdatool的使用，窗函数，matlab生成数字滤波器，通过stm32f103对滤波器进行仿真。其中对fdatool的使用和stm32的使用在其他文章中进行了简述。末文附上学习过程中的参考链接，非常感谢各位大佬的学习分享。
复制链接

扫一扫