几种常用信号平滑去噪的方法(附Matlab代码) |
您所在的位置:网站首页 › 常用实时动态实现的方法 › 几种常用信号平滑去噪的方法(附Matlab代码) |
几种常用信号平滑去噪的方法(附Matlab代码)
1 滑动平均法1.0 移动平均法的方法原理1.1 matlab内自带函数实现移动平均法1.2 利用卷积函数conv()实现移动平均法1.3 利用filter滤波函数实现移动平均法1.4 移动平均的幅频响应1.5 时域和频域的转换关系
2 Savitzky-Golay法2.1 Savitzky-Golay法的方法原理2.2 Savitzky-Golay法的matlab实现2.3 Savitzky-Golay法的幅频响应
3 处理离群值(粗大误差)的方法3.1 中位值法3.2 标准差法和MAD法3.3 Matlab中其它离群值消除方法
4 其它一些FIR滤波器实现光滑去噪4.1 FIR和IIR的区别4.2 利用Matlab构建FIR滤波器
5 IIR滤波器实现光滑去噪6 小波去噪方法
2020年8月更新:增加一个时域和频域的转换关系图,增加相应小节1.5 2021年8月更新:增加第6章小波去噪的内容 信号在实际测量中,难免会混入各种噪声。通常我们希望去除高频的随机噪声,或者是偏离正常测量太大的离群误差,以获得低频的测量数据。下面介绍几种常用的信号平滑去噪的方法。 还是惯例声明:本人没有相关的工程应用经验,只是纯粹对相关算法感兴趣才写此博客。所以如果有错误,欢迎在评论区指正,不胜感激。本文主要关注于算法的实现,对于实际应用等问题本人没有任何经验,所以也不再涉及。 1 滑动平均法本章参考目录 1《数字信号处理》-胡广书 2《Digital Signal Processing - A Practical Guide For Engineers and Scientists》 - Steven W.Smith 1.0 移动平均法的方法原理作为开篇第一个方法,会夹带一些数字信号处理的基本方法,可能会导致篇幅比较啰嗦,之后几章我会尽量挑重点的讲解。 滑动平均法(moving average)也叫做移动平均法、平均法、移动平均值滤波法等等,是一种时间域思想上的信号光滑方法。算法思路为,将该点附近的采样点做算数平均,作为这个点光滑后的值。 matlab有两个函数实现滑动平均法,一个是smoothdata()函数,一个是movmean()函数。 以窗口长度为5为例,smoothdata()函数调用方法为: y = smoothdata( x , 'movmean' , 5 );但是这个smoothdata函数实际上是调用了movmean()函数。所以如果直接使用的话,直接用movmean()会更快。 movmean()函数的调用方法为: y = movmean( x , 5 );下面以一个加噪声的正弦信号为例: %移动平均滤波 clear clc close all N_window = 5;%窗口长度(最好为奇数) t = 0:0.1:10; A = cos(2*pi*0.5*t)+0.3*rand(size(t)); B1 = movmean(A,N_window); figure(1) plot(t,A,t,B1)根据之前的公式,我们可以看到 y ( n ) = 1 / 3 ∗ ( x ( n − 1 ) + x ( n ) + x ( n + 1 ) ) y(n)=1/3*( x(n-1)+x(n)+x(n+1) ) y(n)=1/3∗(x(n−1)+x(n)+x(n+1)) 就相当于一个对x和向量[1/3 1/3 1/3]做卷积。 可以验证: N_window = 5;%窗口长度(最好为奇数) t = 0:0.25:10;%时间 A = cos(2*pi*0.5*t)+0.3*rand(size(t)); B1 = movmean(A,N_window); F_average = 1/N_window*ones(1,N_window);%构建卷积核 B2 = conv(A,F_average,'same');%利用卷积的方法计算 figure(2) plot(t,B1,t,B2)
首先介绍一下Z变换。以向前的滑动平均为例(这里中间值不是n而是n+1,所以相位会移动)。 y ( n ) = 1 / 3 ∗ ( x ( n ) + x ( n + 1 ) + x ( n + 2 ) ) y(n)=1/3*( x(n)+x(n+1)+x(n+2) ) y(n)=1/3∗(x(n)+x(n+1)+x(n+2)) 它的Z变换可以简单的理解为,把x(n+k)替换为z^(-k),即 H ( z ) = 1 3 ( z 0 + z − 1 + z − 2 ) H(z)=\frac{1}{3}(z^{0}+z^{-1}+z^{-2}) H(z)=31(z0+z−1+z−2) 因此对于filter滤波函数,输入的格式为: y = filter(b,a,x)其中b和a的定义为: H ( z ) = b 1 + b 2 ⋅ z − 1 + b 3 ⋅ z − 2 + ⋯ + b n ⋅ z − n + 1 1 + a 2 ⋅ z − 1 + a 3 ⋅ z − 2 + ⋯ + a m ⋅ z − m + 1 H(z)=\frac{b_1+b_2\cdot z^{-1}+b_3\cdot z^{-2}+\cdots +b_{n}\cdot z^{-n+1}}{1+a_2\cdot z^{-1}+a_3\cdot z^{-2}+\cdots +a_{m}\cdot z^{-m+1}} H(z)=1+a2⋅z−1+a3⋅z−2+⋯+am⋅z−m+1b1+b2⋅z−1+b3⋅z−2+⋯+bn⋅z−n+1 其中a1必须为1。所以对应的移动平均滤波可以表示为: y = filter( [1/5,1/5,1/5,1/5,1/5] , [1] , x );它和下面代码的是等价的(在边缘上的处理方式有所不同) y = movmean( x , [4,0] );代表有偏移的滑动平均滤波,4是向后4个点的意思,0是向前0个点的意思。 因为 filter滤波器使用有偏移的向后滤波。滤波后,相位会发生改变。所以通常采用零相位滤波器进行滤波,matlab内的函数为filtfilt()。原理从函数名字上就可以体现出来,就是先正常滤波,之后再将信号从后向前再次滤波。正滤一遍反滤一遍,使得相位偏移等于0。这个之后再IIR滤波器会讲一下。 1.4 移动平均的幅频响应幅频响应可以通过之前1.3得到的H(z)函数来得到,在单位圆上采样,也就是把z替换为e^iw。 以中心窗口为例, y ( n ) = 1 3 ( x ( n − 1 ) + x ( n ) + x ( n + 1 ) ) y(n)=\frac{1}{3} ( x(n-1)+x(n)+x(n+1) ) y(n)=31(x(n−1)+x(n)+x(n+1)) H ( i w ) = 1 3 ( exp ( i w ) 1 + exp ( i w ) 0 + exp ( i w ) − 1 ) H(iw)=\frac{1}{3}(\exp(iw)^{1}+\exp(iw)^{0}+\exp(iw)^{-1}) H(iw)=31(exp(iw)1+exp(iw)0+exp(iw)−1) H(iw)的绝对值就是该滤波方法的幅频响应。以3点滤波为例,matlab代码为: %计算频率响应 N_window = 3; w = linspace(0,pi,400); H_iw = zeros(N_window,length(w)); n=0; for k = -(N_window-1)/2:1:(N_window-1)/2 n = n+1; H_iw(n,:) =1/3* exp(w.*1i).^(-k);%公式(e^iw)^(-k) end H_iw_Sum = abs(sum(H_iw,1));%相加 figure() plot(w/2/pi,H_iw_Sum)由于H变换在单位圆上的特性相当于傅里叶变换,所以上面代码等效于下面计算方法: %计算频率响应 N_window = 3; Y=zeros(400,1); Y(1:N_window) = 1/3;%设置滤波器 Y_F = fft(Y); figure() plot(linspace(0,0.5,200),abs(Y_F(1:200)));matlab也有自带的函数来看频率特性,freqz(),推荐使用这种。 其中,归一化频率等于信号频率除以采样频率f/Fs,采样频率等于时间采样间隔的倒数1/dt。对比不同窗口长度的幅频响应,可以看到: 1平均所采用的点数越多,高频信号的滤波效果越好。 2 3点平均对于1/3频率的信号滤波效果最好,5点平均对1/5和2/5频率的信号滤波效果最好。所以根据这个特性,一方面我们要好好利用,一方面也要避免其影响。 反之,以美国近期的疫情为例,疫情的采样频率为1天一采样,而且显示出很强的7日一周期的特性(病毒也要过周末)。所以,为了消除这个归一化频率为1/7的噪声影响,采样7点的滑动平均滤波。可以看到所有以7天为一变化的信号分量全部被消除掉了。 (下面这个图经常被引用,但是很少有人思考为什么用7天平均方法来平滑数据。) 所以总结一下,移动平均滤波拥有保低频滤高频的特点,而且对于特点频率的滤波具有良好的效果。但是缺点是所有频率分量的信号都会有不同程度衰减。 1.5 时域和频域的转换关系额外补充一部分小内容,可能前面有些概念加入的太突然。很多人可能觉得之前时域上的平均法非常好理解,为什么突然加入幅频特性图,又是Z变换又是fft的。 其实时域上的滤波和频域上的滤波是可以互相转换,且一一对应的。也就是时域上的卷积等于频域上的乘积。 下图为3点移动平均滤波法,时域和频域的转换关系: 虽然前面的 movmean()或者conv()等函数都是用时域实现的信号滤波,但是同样也可以完全在频域上实现。采用ifft(fft(x).*fft(F))实现的滤波效果,和完全时域上的滤波效果是等价的。 下面是展示了窗口长度为3的平滑滤波,从时域上和频域上对信号进行滤波的对比: %实验,检验频域和时域的一致性 %以3点滤波为例 clear clc close all N_window = 3;%窗口长度(最好为奇数) t = 0:0.1:10; A = cos(2*pi*0.3*t)+0.1*cos(2*pi*5*t)+0.2*randn(size(t)); F_average = 1/N_window*ones(1,N_window);%创建滤波器 B2 = conv(A,F_average,'same');%利用卷积的方法计算 figure(1) plot(t,A,'k','linewidth',0.8) %计算原信号的fft A_fft=fft(A); %构建频域上的滤波器 F_average2=zeros(size(t));%长度与x相同,为了后面.*运算 F_average2(1:(N_window-1)/2+1) = 1/N_window; F_average2(end-(N_window-1)/2+1:end) = 1/N_window;%前后设置对称,使得相位不变 F_Fft = fft(F_average2); figure(2) subplot(1,2,1) plot(linspace(0,1,length(t)),abs(A_fft),'linewidth',1); subplot(1,2,2) plot(linspace(0,1,length(t)),abs(F_Fft),'linewidth',1); %进行反逆变换 B3=ifft(A_fft.*F_Fft); figure() plot( t,B2,t,B3 )%对比时域操作和频域操作的差异这也意味着你也可以在频域上操作,实现想要的滤波。比如想要低频通过高频衰减,就把fft后的信号,高频部分强行等于0即可。比如想要消除某个频率的信号(陷波),就令fft后那个信号的频率等于0即可。同理,想要把振幅衰减1/2,就在对应频域上乘以0.5。 2 Savitzky-Golay法本章参考目录: 1《数字信号处理》-胡广书 2 Savitzky-Golay平滑滤波器的最小二乘拟合原理综述 https://wenku.baidu.com/view/b63017eed0d233d4b04e690e.html?fr=search 2.1 Savitzky-Golay法的方法原理Savitzky-Golay法,又叫做平滑滤波器,最著名的就是5点3次滤波器。这是一种基于时间域上的多项式拟合,来消除噪声的方法。 以简单的5点2次构造方法为例,介绍一下基本的求解系数的方法: 首先选取5个点x[-2],x[-1],x[0],x[1],x[2],根据这5个点,构造一条2次抛物线: f ( i ) = a 0 + a 1 ⋅ i + a 2 ⋅ i 2 f(i)=a_{0}+a_{1} \cdot i+a_{2} \cdot i^2 f(i)=a0+a1⋅i+a2⋅i2 这里i=-2,-1,0,1,2。因此,我们需要寻找最优的a0,a1,a2,使得最小二乘拟合最小。最小二乘拟合的函数为: E = ∑ ( f ( i ) − x ( − 2 ) ) 2 = ( f ( − 2 ) − x ( − 2 ) ) 2 + ( f ( − 1 ) − x ( − 1 ) ) 2 + + ( f ( 0 ) − x ( 0 ) ) 2 + ( f ( 1 ) − x ( 1 ) ) 2 + ( f ( 2 ) − x ( 2 ) ) 2 E=\sum{(f(i)-x(-2))^2}\\ =(f(-2)-x(-2))^2+(f(-1)-x(-1))^2+\\ +(f(0)-x(0))^2+(f(1)-x(1))^2+(f(2)-x(2))^2 E=∑(f(i)−x(−2))2=(f(−2)−x(−2))2+(f(−1)−x(−1))2++(f(0)−x(0))2+(f(1)−x(1))2+(f(2)−x(2))2 之后希望最小二乘E最小,我们使得其导数等于0,也就是 ∂ E ∂ a p = 0 \frac{\partial E}{\partial a_p}=0 ∂ap∂E=0 上面式子中,a0、a1和a2三个偏导数,可以得到3个方程。联立,即可求得a0、a1和a2。 对于无相位差的滤波,我们希望窗口是对称的。也就是说,我们希望用5个点,去估计f[0]的值。因此我们需要的只有a0,因为: f ( 0 ) = a 0 + a 1 ⋅ 0 + a 2 ⋅ 0 = a 0 f(0)=a_{0}+a_{1} \cdot 0+a_{2} \cdot 0=a_{0} f(0)=a0+a1⋅0+a2⋅0=a0 这里借助Mathematica的符号运算进行编程求解,程序如下: Clear["Global`*"] M = 5;(*M个点*) P = 2;(*P次*) Mk = (M - 1)/2; f[i_] = Sum[a[k]*i^k, {k, 0, P}] En = Sum[(f[k] - x[k])^2, {k, -Mk, Mk}]; DEn = D[En, {Array[a, P + 1, {0, P}]}]; s = Solve[DEn == 0, Array[a, P + 1, {0, P}]]; s[[1]][[1]]可以得到结果: a[0] -> 1/35 (-3 x[-2] + 12 x[-1] + 17 x[0] + 12 x[1] - 3 x[2])对于经典的5点3次平滑结果,只需要把上面代码P改为3即可,恰好结果相同: a[0] -> 1/35 (-3 x[-2] + 12 x[-1] + 17 x[0] + 12 x[1] - 3 x[2])Savitzky-Golay法可以看做是移动平均法的推广。因为当次数P等于1时,Savitzky-Golay法退化为求解平均值的形态。当然,它也可以视为一种加权平均。其它加权平均比如高斯加权等等,可以在smoothdata函数里的帮助文档看到,就不细说了。 除此之外,也可以直接用matlab中的sgolay()函数,直接得到系数,还是以5点3次为例: M=5;%5点 P=3;%3次 b = sgolay( P , M ); a0 = b((M+1)/2,:) 2.2 Savitzky-Golay法的matlab实现matlab中实现方法和第一章中滑动滤波的方法相似。 利用matlab自带的smoothdata(A,‘sgolay’)函数就可以实现Savitzky-Golay法滤波。但是,该函数只支持N点2次的滤波(2.1阶证明了5点2次和5点3次的a0滤波系数相同,但是其他点数则不相同,这一点一定注意)。 如果想得到不同点数不同次数的公式,参考2.1中的sgolay()函数,可以得到不同的系数。下面代码演示了一个利用matlab自带的5点2阶插值和自己编程实现的5点3次插值,其中自己编程的部分采用了1.2节部分的卷积运算,和2.1节介绍的matlab自带函数sgolay()。 %sgolay滤波 clear clc close all N_window = 5;%窗口长度(最好为奇数) t = 0:0.1:10; A = cos(2*pi*0.5*t)+0.4*rand(size(t)); %matlab自带的5点2次插值 B1 = smoothdata(A,'sgolay' ,N_window); figure(1) plot(t,A,t,B1) %5点3次插值 B2 = smooth_SG_hyh(A,5,3,1); figure(2) plot(t,B1,t,B2) function y_new = smooth_SG_hyh(y,M,P,n) %M点P次插值 %M为窗口长度 %P为拟合阶数 %n为光滑n次 m = length(y); N_window = M;%窗口长度 b = sgolay( P , N_window ); F_SG = b((N_window+1)/2,:);%5点3次插值 y_new = y; for k=1:n y_new2 = wextend(1,'sym',y_new,(N_window-1)/2);%镜像延拓 y_new2 = conv(y_new2,F_SG,'same');%利用卷积的方法计算 y_new = wkeep1(y_new2,m);%中间截断 end end利用之前1.4的内容,可以求出Savitzky-Golay法的幅频响应。 %频率响应对比 [w,A_W_5_3] = A_W_SG_hyh(5,3); [w,A_W_5_1] = A_W_SG_hyh(5,1); figure() plot(w,A_W_5_1,w,A_W_5_3) %频率响应对比 [w,A_W_7_3] = A_W_SG_hyh(7,3); [w,A_W_7_1] = A_W_SG_hyh(7,1); figure() plot(w,A_W_7_1,w,A_W_7_3) ylim([0,1]) function [w_2p,H_iw_Sum] = A_W_SG_hyh(M,P) N_window=M; b = sgolay( P , N_window ); F_SG = b((N_window+1)/2,:);%5点3次插值 w = linspace(0,pi,400); H_iw = zeros(N_window,length(w)); n=0; for k = -(N_window-1)/2:1:(N_window-1)/2 n = n+1; H_iw(n,:) =F_SG(n)* exp(w.*1i).^(-k);%公式(e^iw)^(-k) end H_iw_Sum = abs(sum(H_iw,1));%相加 w_2p = w/2/pi; end
离群值是指在测量值中,出现了一些反常的值,这个反常值与附近的其它正常值差别非常大。 常用的方法有利用 3 σ 3 \sigma 3σ、中位值法、以及基于中位值的MAD方法等等。 本章参考目录: 1 离群值处理方法 https://blog.csdn.net/zjz199303/article/details/79135530 2 常用去除离群值的算法 https://blog.csdn.net/dulingwen/article/details/97006884 3.1 中位值法中位值法,也叫移动中位数法、中值滤波法等。其思想是将窗口内的中位数作为输出结果,如下图所示: matlab中自带的movmedian()函数可以实现中位值法滤波。 clear clc close all %离群值的删除 t = 0:0.05:10; A = sin(t); Ri = randi(length(t),10,1); A(Ri) = A(Ri)*2; %移动中位数,窗口设置为7 B1 = movmedian(A,7); figure(1) plot(t,A,t,B1)下图可以看到中位值法对于离群值处理的结果。 标准差法又叫做 3 σ 3 \sigma 3σ法。目的是规定一个数据波动阈值,当数据超过这个阈值的时候,便认为该数据离群。这个方法阈值的选取方法,采用窗口数据的3倍标准差。 MAD法也是定义了一个阈值,这个阈值叫做中位数绝对偏差MAD。如果超过了3倍的MAD,则认为该数据离群。 两者在Matlab里,可以用filloutliers()函数进行实现。下面代码对比了标准差法和MAD法在消除离群值的效果。 clear clc close all %离群值的删除 t = 0:0.06:10; A = sin(t)+0.1*rand(size(t)); Ri = randi(length(t),4,1); A2 = A; A2(Ri) = A(Ri)*3; figure(3) B2 = filloutliers(A2,'linear','movmean',11);%标准差法 B3 = filloutliers(A2,'linear','movmedian',11);%MAD法 subplot(3,1,1) plot(t,A2) YL = ylim; subplot(3,1,2) plot(t,B2) ylim(YL) subplot(3,1,3) plot(t,B3) ylim(YL)下图可以看到,原本4个离群值,标准差法只找到了1个,但是MAD方法能够做到全部找到,说明MAD方法比标准差法的效果更好,一般更推荐用MAD方法。 matlab自带的函数smoothdata还有两种方法,可以兼顾去噪和去除离群噪声,一个是’rlowess’ 方法,一个是’rloess’ 方法。下面以’rloess’ 方法为例: t = 0:0.06:10; A = sin(t)+0.2*rand(size(t)); Ri = randi(length(t),4,1); A2=A; A2(Ri) = A(Ri)*3; figure(4) B4 = smoothdata(A2,'rloess',8) ; plot(t,A2,t,B4) legend('原函数','光滑')下图可以看到,数据在光滑的同时,顺便也把离群值去掉了。 本章参考: Matlab数字滤波入门 https://zhuanlan.zhihu.com/p/65483011 4.1 FIR和IIR的区别FIR滤波器也叫做 有限冲击响应滤波器。和它相对的是IIR 无限冲击响应滤波器。在matlab里,有一张图可以比较直观。 对于FIR滤波器,由于只有分子b,所以可以用卷积conv()函数去进行滤波。但是对于IIR,就不能用卷积函数去计算。一般认为,IIR函数在相同阶数下,滤波效果要比FIR函数要好,但是有可能出现发散问题。 比如前面的移动平均滤波中的1.3章所介绍的,b=[1/3,1/3,1/3],a=1。就属于一个典型的FIR滤波器。其中,conv(x,b)卷积方法相当于无相位偏移的中心滤波, y ( n ) = 1 3 ( x ( n − 1 ) + x ( n ) + x ( n + 1 ) ) y(n)=\frac{1}{3} ( x(n-1)+x(n)+x(n+1) ) y(n)=31(x(n−1)+x(n)+x(n+1)) 而filter(b,a,x)函数相当于有相位偏移的向后滤波 y ( n ) = 1 3 ( x ( n ) + x ( n + 1 ) + x ( n + 2 ) ) y(n)=\frac{1}{3} ( x(n)+x(n+1)+x(n+2) ) y(n)=31(x(n)+x(n+1)+x(n+2)) ; 事实上,通过下面的等式,我们还可以构建一个等价的IIR滤波器: H ( z ) = 1 3 ( z 0 + z − 1 + z − 2 ) = 1 3 1 − z − 3 1 − z − 1 H(z)=\frac{1}{3}(z^{0}+z^{-1}+z^{-2})=\frac{1}{3}\frac{1-z^{-3}}{1-z^{-1}} H(z)=31(z0+z−1+z−2)=311−z−11−z−3 此时b=[1,0,0,0,-1],a=[1,-1]。 Fs=20; t = 0:1/Fs:10; A = cos(2*pi*0.3*t)+1*sin(2*pi*1.2*t)+0.8*rand(size(t)); b=[1/3,1/3,1/3];a=1; B1 = filter(b,a,A); b=1/3*[1,0,0,0,-1];a=[1,-1]; B2 = filter(b,a,A); plot(t,B1,t,B2)对比这两个滤波器,可以看到滤波效果是基本等价的。 回到FIR滤波器上,之前提到的时域角度提出的均值滤波器和Savitzky-Golay滤波器,都属于FIR滤波器。他们在频域上的表现并不能随意的控制。 如果我们希望构建一个滤波器,让它在频域上满足自己的要求,则可以用fir1()函数去构建。我这里只演示一个简单地,更专业的就不在这里献丑了。 比如,我的采样频率为10hz,我想提取的信号在1hz附近,而噪声则大于等于2hz。因此,需要构建一个低频通过,高频减弱的滤波器,使得在1.5hz以下的信号能够保留,1.5hz以上的信号删除。 于是用下面代码去构建,滤波器长度为31,保留0-0.15之内的信号。(注1:matlab的归一化信号不是0-0.5,而是0-1,所以对应的0.15需要乘以2倍,变成0.3)(注2:最小频率不能是0,必须大于0,所以用0.001代替)(注3:滤波器长度为31是为了演示,实际上不需要这么大的长度,滤波器长度大了的话计算量也会增大) b = fir1(31,[0.001,0.3]);构建的滤波器形状和频率特性为: 频率特性满足预期的设计。滤波效果如下: 比较著名有巴特沃斯滤波器butter()和切比雪夫滤波器cheby1()。这里感觉还有地方没看懂,就不详细写了。matlab里可以通过designfilt()来自定义的构建滤波器。 以butter滤波器为例,采用零相位滤波,和第4.2节一样,希望保留0-0.15之内的信号。那么构建滤波器格式为: [b,a] = butter(6,0.15*2,'low' );b为分子,a为分母。6为6阶,代表滤波器的长度为6+1个点。'low’是低通滤波器的意思,0.15*2代表保留0-0.15之内的信号。得到的滤波器幅频响应曲线如下,可以看到能够满足要求。 平滑后的曲线效果如下图。 小波去噪最近刚刚接触,这里也更新补充一下。后续具体原理和算法,我尽快补充上(又挖坑)。和传统去噪相比,其优势在于可以更好的提取出信号的主要部分,受噪声影响较小。 比如下面这个图,传统的移动平均去噪(也包括前面介绍的各种滤波器),由于在频域上相当于全局频域上的滤波,所以会出现:高频信号和噪声一块被滤掉了,中间频域去噪效果比较好,低频受噪声影响严重滤波不足 的现象。这就是全局滤波的缺点,按下葫芦浮起瓢,不能做到高频低频同时兼顾。 对于非单一频率非周期变化的信号来说,小波去噪的效果明显更好。 matlab有自带的小波去噪函数wden()和wdenoise()。代码如下: clear clc load noisdopp y1 = wdenoise(noisdopp); y2 = smoothdata( noisdopp , 'movmean' , 5 ); N=length(y1); figure() subplot(2,1,1) hold on plot(1:N,noisdopp) plot(1:N,y1,'linewidth',1) hold off xlim([0,1000]) legend('原信号','小波去噪','Location','southeast') subplot(2,1,2) hold on plot(1:N,noisdopp) plot(1:N,y2,'linewidth',1) hold off xlim([0,1000]) legend('原信号','移动平均','Location','southeast') |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |