IIR工频陷波器问题

天陨石 · 发表于 2014-8-19 14:53:31

最近在弄一个项目，刚好要用IIR工频陷波器，看到网上这方面的资料不是很多，自己完成之后就在博客中总结了一下，以前老是在原子哥的网站上拿东西，现在有时间写了下，第一次发帖，希望和大家分享，废话不多说了。

工频干扰是由电力系统引起的一种干扰，中国大陆地区由50hz及其谐波构成的一种干扰，幅值约为ECG（心电信号的50%）。

工频干扰可以用陷波器进行处理，过滤50hz信号。常用的工频陷波器主要有IIR和FIR两种，其中FIR具有良好的线性相位，但是在同等滤波效果的情况下，IIR的阶数要比FIR少很多，一个两阶的IIR滤波器的效果FIR要付出100多阶的代价，阶数大意味着运算量大，对于一个MCU单片机来说这是得不偿失的，所以采用IIR滤波器来实现工频滤波。

IIR滤波器的设计方法有脉冲响应不变法和双线性变换法，如果你已经熟悉IIR滤波器的设计原理可以直接借用Matlab的Fdatool滤波器设计工具直接生成相关系数，这样省去了中间设计环节，再通过IIR滤波器的直接II型实现方式翻译成C语言格式并应用于STM32平台（重要：fs为采样频率，Fnotch为陷波频率）

按Designer键即可生成IIR滤波系数

导出滤波系数得到转移函数

这样再通过直接II型实现框图，翻译成C语言代码并烧录到主控MCU中运行

当然你也可以选择直接I型实现，只不过直接II型比直接I型少占一半的内存空间，所以选择直接II型实现。

[mw_shl_code=c,true] x0=ADC_ConvertedValueLocal; //输入ADC采集到的信号 w0[0]=IIR_50Notch_A[0]*x0-IIR_50Notch_A[1]*w0[1]-IIR_50Notch_A[2]*w0[2]; y0=IIR_50Notch_B[0]*w0[0]+IIR_50Notch_B[1]*w0[1]+IIR_50Notch_B[2]*w0[2]; w0[2]=w0[1]; w0[1]=w0[0]; [/mw_shl_code]

前提是你的输入信号的时域采样率要和滤波器的频域采样率要保持一致，这样就涉及定时器触发ADC采样，用STM32的定时器来配置采样率，

你可以先用Matlab平台来验证你滤波器设计是否正确

附带Matlab平台仿真实现：

[mw_shl_code=c,true]clear all fs=250; %数字滤波器的采样频率fs=1000hz f=50; %50hz的正弦信号 t=0:1/fs:4; %f为时间间隔，ADC采样频率250hz s=sin(2*pi*f*t); IIR_B=[0.90239774423695518,-0.55771247730967288,0.90239774423695518]; IIR_A= [1,-0.55771247730967288,0.80479548847391036]; w01=0; w02=0; w03=0; for i=1:1000 w01=s(i)-IIR_A(2)*w02-IIR_A(3)*w03; y0(i)=IIR_B(1)*w01+IIR_B(2)*w02+IIR_B(3)*w03; w03=w02; w02=w01; end figure subplot(211) axis([1,1000,-1,1]); hold on plot(s); title('直接II型实现—原始信号时域'); subplot(212) plot(y0,'r'); title('直接II型实现—信号滤波之后'); [/mw_shl_code]

附带一STM32实现IIR滤波器的源代码:

http://pan.baidu.com/s/1dDcRCKp 提取密码为：muye

zuozhongkai · 发表于 2014-8-19 15:10:52

牛！！！！！！！！学习了

正点原子 · 发表于 2014-8-19 23:06:22

不错,谢谢分享.

cornrn · 发表于 2014-12-9 16:58:37

这才是有料的帖子，谢谢分享

cornrn · 发表于 2014-12-10 08:36:38

我用的MATLAB2012b怎么没见到陷波滤波器啊？

cornrn · 发表于 2014-12-10 10:33:30

还有就是这个转移函数分母上1怎么被你省了？

cornrn · 发表于 2014-12-10 10:33:55

回复【2楼】zuozhongkai:
---------------------------------
拉住你们关注帖子

cornrn · 发表于 2014-12-10 10:34:09

回复【3楼】正点原子:
---------------------------------
拉住你们关注帖子。。。

Di诺 · 发表于 2014-12-12 11:14:03

看不懂。

cornrn · 发表于 2014-12-22 18:29:31

回复【9楼】Di诺:
---------------------------------
做出来了，感谢这帖主

Huangogo · 发表于 2015-4-15 16:47:49

回复【10楼】cornrn:
---------------------------
你好，请问你是怎么传递系数的呀？我也在用IIR，可是matlab生成的系数不对

q1092840790 · 发表于 2015-4-29 21:45:28

楼主，你用matlab设计的采样频率为250Hz，而在代码中RCC_ADCCLKConfig(RCC_PCLK2_Div6); ADC_RegularChannelConfig(ADC1, ADC_Channel_10, 1, ADC_SampleTime_1Cycles5);即ADC时钟为12MHz，转换周期为14，因此采样频率应该是6/7MHz才对，怎么会是250Hz？
          对于结构系数const float IIR_50Notch_B[3] = {
                                             0.9023977442,
                                    -0.5577124773,
                                     0.9023977442
                                         };
const float应该只有6位小数位有效...这明显不对额？

q1092840790 · 发表于 2015-4-29 21:53:49

回复【10楼】cornrn:
---------------------------------
朋友，能不能帮忙看下楼上我的问题？

q1092840790 · 发表于 2015-4-30 12:06:39

以上两问已理解。。。不知为什么要在代码里开启PWM？。。。

q1092840790 · 发表于 2015-4-30 13:39:30

而且这个TIM2 CH2对应的GPIO口PA1都没有进行初始化

cornrn · 发表于 2015-4-30 19:51:46

回复【12楼】q1092840790:
---------------------------------
你确定看到我程序了？我用SYSTICK中断，确定采样时间的，4ms进systick开启ADC一次扫描就可以确保250HZ采样率了，再说谁那么傻一个点持续那么长时间采，那么那个值还算真值吗，理论上采集时间越短才能确保是那个时间节点上的当前值，不然就是T1~T2时间上的均值

q1092840790 · 发表于 2015-5-1 18:28:14

回复【16楼】cornrn:
---------------------------------
- -我说的是贴主1楼的那个程序。。。。采样时间那个是我的问题，没理解好，不好意思....
我用time中断触发ADC采样了，采样功能已经实现了。发现250Hz采样频率根本行不通，对于50Hz工频而言，一周期就采5个点，并没有什么用。。。采样出来的不经处理的波形跟矩形波差不多。用这个帖子里的处理方法，实际上效果基本没有
我现在在做IIR滤波......50Hz 工频陷波器，50Hz 时信号衰减大于等于60dB，48Hz 时信号衰减≤0.1dB，52Hz 时信号衰减≤0.1dB；用直接I型高阶IIR没弄出来。。。。麻烦问你你是怎么实现的？

天陨石 · 发表于 2015-5-2 10:42:47

回复【17楼】q1092840790:
---------------------------------
1. 采样频率50hz,一个周期采5个点没什么用
答：可以加大采样周期，具体操作如下：matlab fdatool设计滤波器时，设定你想设置的采样频率。对应stm32 时钟触发ADC采样的频率需要和滤波器一致。（ps: 1000hz的采样频率为佳）
2. 50hz时信号衰减大于等于60db， 48hz衰减信号小于等于0.1, 52hz信号衰减小于等于0.1db
答：可以不用IIR陷波器，用fdatool设计IIR带阻滤波器，上通带限为52hz，下通带限为48hz, 阻带衰减为60dB,通带衰减为0dB。选择滤波器设计为最小阶数。
3. 关于陷波器和带阻滤波器
答：陷波器其实质是阻带趋于点的带阻滤波器
4. 关于IIR直接I型和直接II型
答：尽量采用直接II型。直接II型省了一半的内存量。当然由于IIR阶数不高，所以两者在内存消耗上没有什么区别。
5. 关于fdatool设计的滤波器怎么用C语言实现
答：5.1 首先通过fdatool得到你滤波器的抽头系数。
5.2 通过matlab平台，仿真检测滤波器实现方式的正确性。（通过编写.m脚本文件）一般参数设置正确，抽头系数不会有问题。问题一般在实现方式
5.3 先产生一个50hz的正弦波（注意时域采样频率需要和滤波器采样频率一致），然后通过你的滤波器（可以用for循环实现，参照帖子上面matalb滤波器实现部分）,查看输出，纵坐标单位设置为dB
5.4 产生一个48hz的正弦波，然后通过你的滤波器
5.5 产生一个52hz的正弦波，然后通过你的滤波器
5.6 matlab仿真没问题之后，可以翻译成stm32 c语言
5.7 如果存在问题，这时候肯定不会是滤波器问题，请查看定时器触发ADC采样等部分代码

q1092840790 · 发表于 2015-5-2 11:10:54

回复【18楼】天陨石:
---------------------------------
你好！感谢你的回复。我将采样频率设为了5000Hz，即分频系数为36-1，重载数为400-1.由于采用的是带阻滤波器，因此不可能一组数据就能算完，需要用到数组。
float iir_fliter(float in,const float IIR_B[][3],const float IIR_A[][3],int IIR_NSEC)
{
int i;
x[0][0]=in;
for(i=0;i<IIR_NSEC;i++)
{
y[0]=x[0]*IIR_B[0] + x[1]*IIR_B[1] + x[2]*IIR_B[2] - y[1]*IIR_A[1] - y[2]*IIR_A[2];
y[0]/=IIR_A[0];

y[2]=y[1];y[1]=y[0];
x[2]=x[1];x[1]=x[0];

x[i+1][0]=y[0];
}

if(x[IIR_NSEC][0]>4096) x[IIR_NSEC][0]=4096;
if(x[IIR_NSEC][0]<0) x[IIR_NSEC][0]=0;

return (x[IIR_NSEC][0]);
}
我编写了上面这个函数，实现了带阻的功能。但50Hz滤波只是我需要完成的一部分。
还需要做10Hz低通滤波，18Hz衰减大于40dB；1Hz高通滤波，0.5Hz衰减>40dB.
当我用fdatool设计出系数时，继续使用上面这套算法，发现只有带阻能够实现，高通和低通出来的波形均不对，是否上面这套算法不适用于低通和高通？我使用的是IIR直接II型

q1092840790 · 发表于 2015-5-2 12:31:55

回复【18楼】天陨石:
---------------------------------
高通未能实现原因是：高通滤波器会滤去直流成分，导致输出波形含有负分量，我在输出的的基础上加了一直流偏置，实现了高通滤波。。。。
目前低通滤波一直无法实现，示波器在直流耦合的情况下，输出都为0，无论频率如何改变

天陨石 · 发表于 2015-5-2 14:24:28

回复【20楼】q1092840790:
---------------------------------
1HZ高通滤波就是抵制低频成分，直流是0频分量属于被抑制范围

q1092840790 · 发表于 2015-5-2 14:26:01

回复【21楼】天陨石:
---------------------------------
高通我做出来。。。。低通无论如何修改fdatool设计，出来都是0 T.T.求助

天陨石 · 发表于 2015-5-2 14:30:42

回复【22楼】q1092840790:
---------------------------------
1. 看看低通参数是否正确
2. 若参数没有问题，通过matlab仿真，产生一个高频正弦波，通过滤波器。看看是否被抑制
3. 若高频能被抑制，产生一个低频正弦波（通带范围内），通过滤波器。看看是否不被抑制

q1092840790 · 发表于 2015-5-2 14:53:32

回复【23楼】天陨石:
---------------------------------
实际上我并不会使用matlab 。。。我只是在短时间突击内学会了fdatool的简单应用。能不能麻烦您帮我验证下？

q1092840790 · 发表于 2015-5-2 15:57:33

回复【23楼】天陨石:
---------------------------------
刚突击了下。。。。生成Objiect Hd
>> % 1. 产生信号（频率为10Hz和100Hz的正弦波叠加）
Fs=5000; % 采样频率5000Hz
t=0:1/Fs:1;
s10=sin(20*pi*t); % 产生10Hz正弦波
s100=sin(200*pi*t); % 产生100Hz正弦波
s=s10+s100; % 信号叠加
figure(1); % 画图
subplot(2,1,1);plot(s);grid;
title('原始信号');
% 2. FFT分析信号频谱
len = 512;
y=fft(s,len); % 对信号做len点FFT变换
f=Fs*(0:len/2 - 1)/len;
subplot(2,1,2);plot(f,abs(y(1:len/2)));grid;
title('原始信号频谱')
xlabel('Hz');ylabel('幅值');
Bf=filter(Hd1,s); % 进行低通滤波
By=fft(Bf,len); % 对信号f1做len点FFT变换
figure(2); % 画图
subplot(2,1,1);plot(t,s10,'blue',t,Bf,'red');grid;
legend('10Hz原始信号','巴特沃斯滤波器滤波后');
subplot(2,1,2);plot(f,abs(By(1:len/2)));grid;
title('巴特沃斯低通滤波后信号频谱');
xlabel('Hz');ylabel('幅值');

滤波出来的波形是对的。。。。但是运用阻带和高通滤波所用的函数处理参数时，出来的是错误的。。。。全为0

company · 发表于 2015-7-31 16:27:16

楼主分享的代码怎么没有了，麻烦给一份吧谢谢，最近的项目也要陷波，但是实际焊接的确不理想，想利用数字滤波。

天陨石 · 发表于 2015-8-1 21:45:09

回复【26楼】company:
---------------------------------
附件和云盘链接文件一样

FantaSy_ · 发表于 2015-8-1 23:49:14

牛逼啊~~~~

小卿 · 发表于 2016-7-27 18:34:47

楼主，这个50HZ能完全滤掉吗，我现在在用FIR做，都只能滤掉幅值的一半左右，效果还是不行

天陨石 · 发表于 2016-7-28 11:08:13

小卿发表于 2016-7-27 18:34
楼主，这个50HZ能完全滤掉吗，我现在在用FIR做，都只能滤掉幅值的一半左右，效果还是不行

不可能完全滤掉的，只能说让噪声衰减多少dB。
至于FIR的抑制幅值得能力，我在帖子中已经明确提出，一个三阶IIR滤波器的抑制幅值能力相当于1000多阶的FIR滤波器

aaaaalaode · 发表于 2016-7-28 11:24:36

现在还不知道 matlab能做什么只知道那是大神玩的。

小卿 · 发表于 2016-7-28 16:20:33

天陨石发表于 2016-7-28 11:08
不可能完全滤掉的，只能说让噪声衰减多少dB。
至于FIR的抑制幅值得能力，我在帖子中已经明确提出，一个 ...

你好，我也试过把你的滤波程序去试了下，可是发现为什么输入进去250个点，可实际出来的点数却只有100多个点的，其余的点是怎么处理掉了

天陨石 · 发表于 2016-7-30 11:31:38

小卿发表于 2016-7-28 16:20
你好，我也试过把你的滤波程序去试了下，可是发现为什么输入进去250个点，可实际出来的点数却只有100多个 ...

你没有说明怎么使用上面的代码，其次数据的采样率有没有和设计的滤波器采样率是否一致，把代码和思路贴出来。你这样直接抛一个问题，我不知道怎么回答

小卿 · 发表于 2016-7-30 12:52:35

天陨石发表于 2016-7-30 11:31
你没有说明怎么使用上面的代码，其次数据的采样率有没有和设计的滤波器采样率是否一致，把代码和思路贴出 ...

我这边找到问题了，因为我这边手机软件程序支持不了数量很大的处理所有才会有那个问题出现，我想问下，为什么我用你的可以完全滤掉采样率为250，50HZ的信号，但为什么我自己也用matlab把采样率改为1200，我发现改了之后却不能完全滤掉50HZ的信号的（也用过Matlabb仿正过没问题的，可是用到你上面的程序就出现传输到蓝牙就不行了），我想问下这个是什么原因，还有能不能帮我设计一个1200采样率，50hz的陷波参数呀，谢谢了

天陨石 · 发表于 2016-8-16 00:53:44

小卿发表于 2016-7-30 12:52
我这边找到问题了，因为我这边手机软件程序支持不了数量很大的处理所有才会有那个问题出现，我想问下，为 ...

帖子中明确提出，只有滤波器的采样频率和实际信号采样频率一致的时候，滤波器才会发挥效用。设计滤波器用fdatool设计，帖子中也明确给出了设计滤波器的流程

大鸭梨 · 发表于 2016-9-2 18:32:58

楼主请教个问题啊，我就是想做个IIR低通的滤波器，系数是通过File\Export\.....导出的，因为找到这个系数与方程的对应关系，我弄了个50阶的低通感觉代码没什么问题，前面几阶可以输出值，有一定的滤波效果，但阶数高了之后输出就全0了，希望楼主能解惑。
代码有点长。。估计你也不想看。。

大鸭梨 · 发表于 2016-9-2 18:35:26

q1092840790 发表于 2015-5-2 11:10
回复【18楼】天陨石:
---------------------------------
你好！感谢你的回复。我将采样频率设为了5000Hz， ...

可以发一份IIR滤波器的C程序吗？最近在学习编程，遇到问题搞不通。。。。。

大鸭梨 · 发表于 2016-9-2 18:44:54

为什么有的帖子上说用I型，有的用II型？我觉得它们的系数是一样的啊？

LuXing · 发表于 2016-12-9 11:42:34

mark一下

jiney555 · 发表于 2017-1-13 15:29:21

各位dax，滤波后为什么在每帧数据都有一个振荡波形，怎么去掉这个振荡波形，楼主仿真里也有这个波形。

whj467467274672 · 发表于 2018-9-5 09:43:04

666，好资料，学习了，必须关注

whj467467274672 · 发表于 2018-9-11 10:00:03

cornrn 发表于 2014-12-10 10:33
还有就是这个转移函数分母上1怎么被你省了？  

请问一下，这个函数怎么出来的，还是手动推算的

edmund1234 · 发表于 2018-9-11 10:29:38

大鸭梨发表于 2016-9-2 18:32
楼主请教个问题啊，我就是想做个IIR低通的滤波器，系数是通过File\Export\.....导出的，因为找到这个系数与 ...

没写错吧， 50阶？呵呵

天陨石 · 发表于 2018-9-25 10:38:01

whj467467274672 发表于 2018-9-11 10:00
请问一下，这个函数怎么出来的，还是手动推算的

用matlab的工具箱设计的，不用手算。

天陨石 · 发表于 2018-9-25 10:39:01

小卿发表于 2016-7-30 12:52
我这边找到问题了，因为我这边手机软件程序支持不了数量很大的处理所有才会有那个问题出现，我想问下，为 ...

滤波器的采样频率要和数据的采样频率对上，才可以

天陨石 · 发表于 2018-9-25 10:40:10

edmund1234 发表于 2018-9-11 10:29
没写错吧， 50阶？呵呵

阶数太高，过滤后的信号会严重失真

1649921925 · 发表于 2020-10-12 11:08:51

请问IIR的从语言程序怎么编？

帐号		自动登录	找回密码
密码			立即注册

IIR工频陷波器 问题

原子哥极力推荐 /2

IIR工频陷波器问题