第六章 实验数据的平滑滤波

40
60
80
100
120
>> y=ph(x,2,1); >> plot(y)1.5
1
0.5
0
-0.5
-1
-1.5
0
20
40
60
80
100
120
>> y=ph(x,2,2); >> plot(y)1.5
1
0.5
0
-0.5
-1
-1.5
0
20
40
60
80
100
120
min(TA - Y)T (TA - Y)
1 2 4 8 1 1 1 1 T 1 0 0 0 1 1 1 1 1 2 4 8
>> T=[ones(1,5);-2:2;(-2:2).^2;(-2:2).^3]'; >> T*inv(T'*T)*T' ans = 69/70 2/35 -3/35 2/35 -1/70 2/35 27/35 12/35 -8/35 2/35 -3/35 12/35 17/35 12/35 -3/35 2/35 -8/35 12/35 27/35 2/35 -1/70 2/35 -3/35 2/35 69/70
>> t=0:0.1:10; >> x=sin(t)+0.1*(2*rand(1,length(t))-1); >> plot(x)
1.5
1
0.5
0
-0.5
-1
-1.5
0
20
40
60
80
100
120
>> y=ph(x,1,1); 1.5 >> plot(y)
1
0.5
0
-0.5
-1
-1.5
0
20
七点:
>> T=[ones(1,7);-3:3;(-3:3).^2;(-3:3).^3]'; >> T*inv(T'*T)*T' ans = 13/14 4/21 -2/21 -2/21 1/42 4/21 19/42 8/21 1/7 -2/21 -2/21 8/21 19/42 2/7 1/21 -2/21 1/7 2/7 1/3 2/7 1/42 -2/21 1/21 2/7 19/42 2/21 -1/6 -2/21 1/7 8/21
对整个数据:
(5 y1 2 y2 y3 ) / 6 y1 yi ( yi 1 yi yi 1 ) / 3 (i 2,3,..., N 1) y ( y 2 y 5 y ) / 6 N 2 N 1 N N
五点线形平滑(n=2) >> T=[ones(1,5);-2:2]' >> T*inv(T'*T)*T' ans = 3/5 2/5 2/5 3/10 1/5 1/5 0 1/10 -1/5 0
6.3 二次加权移动平均
令: 最小二乘准则:
n t n
2 yit A0 At A t 1 2
min [( A0 A1t A2t 2 ) yi t ]2
五点二次平滑(n=2)
min [( A0 A1t A2t 2 ) yi t ]2
t 2 2
>> A=[1 2 3;4 5 6] A= 1 2 3 4 5 6 >> cumprod(A) ans = 1 2 3 4 10 18 >> cumprod(A,2) ans = 1 2 6 4 20 120
function y = ph(x,n,m) % 数据平滑,2n+1点,m次平滑 if 2*n+1<m+1 error('2n+1应不小于m+1'); end N = length(x); t = (-n:n)'; T1 = ones(2*n+1,1); T2 = repmat(t,1,m); T = cumprod([T1 T2],2); M = T*inv(T'*T)*T';
生成有两个频率的信号 >> t=0:0.1:100; >> y1=sin(2*pi*t)+2*sin(6*pi*t); >> f1=fft(y1); >> w=(1:length(t))/length(t); >> plot(2*w,abs(f1))
900 800 700 600 500 400 300 200 100 0
数据滤波问题:
我们实际获得的各种实验数据或信号中总是存 在各种各样的噪声. 对测量数据或信号进行处 理的过程称为滤波. 滤波分为: 1. 频域滤波 2. 时域滤波 数据平滑是时域滤波的一种.
某物理量y是时间t的函数:
y y(t )
对其进行采样: yi y(ti )
一般是等时间间隔: yi y(iT ) yi 随时间在变,故称为y的时域特征, yi 的傅立叶 变换称为y的频域特征. 由时域特征考虑的滤波 处理叫时域滤波,由频域特征考虑的滤波处理叫 频域滤波.
1/5 1/5 1/5 1/5 1/5
0 1/10 1/5 3/10 2/5
-1/5 0 1/5 2/5 3/5
七点线形平滑(n=3)
>> T=[ones(1,7);-3:3]‘ >> T*inv(T'*T)*T'
ans = 13/28 5/14 1/4 1/7 1/28 -1/14 -5/28
5/14 2/7 3/14 1/7 1/14 0 -1/14 1/4 3/14 5/28 1/7 3/28 1/14 1/28 1/7 1/7 1/7 1/7 1/7 1/7 1/7 1/28 1/14 3/28 1/7 5/28 3/14 1/4 -1/14 0 1/14 1/7 3/14 2/7 5/14 -5/28 -1/14 1/28 1/7 1/4 5/14 13/28
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
设计频域滤波器 >> [h,err,res]=remez(40,[0 0.4 0.48 1],[1 1 0 0]); >> plot(res.fgrid,abs(res.H))
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0
0.1
0.2
0.3
0.4
用它取代yi. n=1: yi=100,152,198,249,318,349,403,452,497,550 yi=150,200,255,305,357,401,451,500
550 500 450 400 350 300 250 200 150 100
1
2
3
4
5
6
7
8
9
10
6.1.2 加权移动平均 m y A At ... A t 令: i t 0 1 m t n, (n 1),..., n 1, n 其中: A0 , A1 ,..., Am 用最二乘法求系数: n 即: 2
>> T=[ones(1,3);-1:1]' T= 1 -1 1 0 1 1 >> T*inv(T'*T)*T' ans = 5/6 1/3 -1/6 1/3 1/3 1/3 -1/6 1/3 5/6
即: yi1 (5 yi 1 2 yi yi 1 ) / 6
yi ( yi 1 yi yi 1 ) / 3 y ( y 2 y 5 y ) / 6 i 1 i i 1 i 1
>> K=medfilt2(J); >> imshow(K)
6.1 实验数据的移动平均
6.1.1 单纯移动平均 采集的N个数据: y1 , y2 ,..., yi ,..., yN 对yi前后对称取2n+1个数据,求其平均值:
n 1 yi yi k 2n 1 k n
i n 1, n 2,..., N n
T T min( Y Y ) ( Y Y ) min( TA Y ) (TA - Y) 矩阵形式: T 1 T A ( T T ) TY 解:
得到:
Y T(TT T)1 TT Y
>> T=[ones(1,5);-2:2;(-2:2).^2]' T= 1 -2 4 1 -1 1 1 0 0 1 1 1 1 2 4
for i=1:n y(i) = M(i,:)*x(1:2*n+1)'; y(N-n+i) = M(n+1+i,:)*x(N-2*n:N)'; end X = []; for i=1:2*n+1 X = [X;x(i:N-2*n-1+i)]; end y(n+1:N-n) = M(n+1,:)*X;
0.5
0.6
0.7
0.8
0.9
1
滤波: >> y2=filter(h,1,y1); >> f2=fft(y2); >> plot(2*w,abs(f2))
900 800 700 600 500 400 300 200 100 0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
图形处理:中值滤波. >> I=imread('eight.tif'); >> J=imnoise(I,'salt',0.02); >> imshow(J)
-1/7 1/14 3/14 2/7 2/7 3/14 1/14
-1/14 0 1/14 1/7 3/14 2/7 5/14
5/42 -1/14 -1/7 -2/21 1/14 5/14 16/21
6.4 三次加权移动平滑
令: 矩阵形式: 五点:
2 3 yit A0 At A t A t 1 2 3
令: 1 2 4 1 1 1 A0 A A T 1 0 0 1 A2 1 1 1 1 2 4
yi 2 yi 2 y y i 1 i 1 Y yi TA Y yi yi 1 yi1 y y i2 i2
min ( yit yi t )
t n
从而求得 yi
yit A0 At 1 令: 最小二乘准则:
n
6.2 线形加权移动平滑
min [( A0 A1t ) yi t ]2
t n
三点线形平滑(n=1):
2 min [( A0 At ) y ] 1 i t t 1 1
>> T*inv(T'*T)*T' ans = 31/35 9/35 9/35 13/35 -3/35 12/35 -1/7 6/35 3/35 -1/7
-3/35 12/35 17/35 12/35 -3/35
-1/7 6/35 12/35 13/35 9/35
3/35 -1/7 -3/35 9/35 31/35
Leabharlann Baidu
2/21 -1/6 -2/21 1/7 8/21 19/42
-1/21 2/21 1/42 -2/21 -2/21 4/21
-1/21
2/21
1/42
-2/21
-2/21
4/21
13/14
repmat(), 由向量生成矩阵. >> x=[1 2 3]; >> repmat(x,2,1) ans = 1 2 3 1 2 3 cumprod(), 生成连积矩阵.
矩阵形式:
令:
1 1 T 1 0 1 1
A0 A A1
yi1 TA Y y i yi1
yi 1 Y y i yi 1
T T 矩阵形式: min(Y - Y) (Y - Y) min(TA - Y) (TA - Y) A (TT T)1 TT Y 解: Y T(TT T)1 TT Y 得到:
七点二次平滑(n=3)
>> T=[ones(1,7);-3:3;(-3:3).^2]'; >> T*inv(T'*T)*T' ans = 16/21 5/14 1/14 -2/21 5/14 2/7 3/14 1/7 1/14 3/14 2/7 2/7 -2/21 1/7 2/7 1/3 -1/7 1/14 3/14 2/7 -1/14 0 1/14 1/7 5/42 -1/14 -1/7 -2/21